## this script creates tables 1&2 and the graphics from the mixed multinomial logit models ##
## note that these graphics will only populate properly if the accompanying Stata ##
## script has already been run and output the proper files... that takes a while ##
## libraries ##

	library(foreign)
	library(mnormt)
	library(klaR) 
	library(plotrix)

## don't forget to set your directory! ##
	setwd('')

## first, the descriptives in table 1 ##

	senators <- read.table('senatorData105to109.txt', sep = '|', header = TRUE)
	summary(senators$unity[senators$party == 'D'])
	summary(senators$unity[senators$party == 'R'])
	sd(senators$unity[senators$party == 'D'])
	sd(senators$unity[senators$party == 'R'])

## read in the data ##
## media-voting unity relationship ##

	articles <- read.table('article info.csv', sep = ',', header = TRUE)
	articles$senator <- as.character(articles$senator)
	unity <- read.dta('senatorUnity.dta')
	unity$senator <- as.character(unity$senator)

## constrain to those senators serving our 10 year sample period ##
	articles <- articles[articles$veteran == 1, ]

## add in icpsr ids and unity scores ##
	articles$unity <- c()
	articles$icpsr <- c()
	for(i in 1:length(articles$senator)){
		articles$unity[i] <- unity$pastUnityMean[unity$senator == articles$senator[i]]
		articles$icpsr[i] <- unity$icpsr[unity$senator == articles$senator[i]]
	}
	
## create focal variable: number of maverick-mentions ##
## articles over total number of articles ##

	articles$mentions <- articles$mav2/articles$all
	
## estimate the model and take a posterior draw ##
	set.seed(0624)
	model <- glm(articles$mentions ~ articles$unity)
	summary(model)
	coef <- rmnorm(n = 10, coef(model), vcov(model))

## and now iterate to reproduce table a5.1 in the appendix ##

	aic <- c()
	for(i in 1:5000){
		data <- articles[round(runif(51, 0.5, 51.49)), ]
		model <- glm(data$mentions ~ data$unity)
		new <- rmnorm(n = 10, coef(model), vcov(model))
		coef <- rbind(coef, new)
		aic[i] <- AIC(model)
	}

	summary(aic)
	sd(aic)
	summary(coef)
	sd(coef[, 1])
	sd(coef[, 2])

## now plot the relationship ##

pdf(file = 'figure1.pdf', width = 10, height = 7)
	plot(1, 1, type = 'n',
		axes = FALSE,
		xlim = c(0.75, 1),
		ylim = c(0.00, 0.20),
		ylab = 'Proportion of Maverick Articles',
		xlab = 'Partisan Unity Score',
		main = 'Effect of Partisan Unity on\nMaverick Coverage in the Media')
	axis(1)
	axis(2)
	points(articles$unity, articles$mentions, pch = 16, col = rgb(1, 0, 0, .2), cex = 1.5)
	text(articles$unity[articles$senator == 'MCCAIN (R AZ)']+0.020, articles$mentions[articles$senator == 'MCCAIN (R AZ)'], 'MCCAIN (R AZ)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'SNOWE (R ME)']+0.020, articles$mentions[articles$senator == 'SNOWE (R ME)'], 'SNOWE (R ME)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'STEVENS (R AK)']+0.020, articles$mentions[articles$senator == 'STEVENS (R AK)'], 'STEVENS (R AK)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'MIKULSKI (D MD)']+0.020, articles$mentions[articles$senator == 'MIKULSKI (D MD)'], 'MIKULSKI (D MD)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'LIEBERMAN (D CT)']-0.022, articles$mentions[articles$senator == 'LIEBERMAN (D CT)'], 'LIEBERMAN (D CT)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'SANTORUM (R PA)']+0.022, articles$mentions[articles$senator == 'SANTORUM (R PA)'], 'SANTORUM (R PA)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'BIDEN (D DE)']+0.015, articles$mentions[articles$senator == 'BIDEN (D DE)']-0.005, 'BIDEN (D DE)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'SPECTER (R PA)']+0.020, articles$mentions[articles$senator == 'SPECTER (R PA)'], 'SPECTER (R PA)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'DOMENICI (R NM)']+0.021, articles$mentions[articles$senator == 'DOMENICI (R NM)'], 'DOMENICI (R NM)', col = rgb(0, 0, 1, .5), cex = .8)
	text(articles$unity[articles$senator == 'FEINGOLD (D WI)']+0.020, articles$mentions[articles$senator == 'FEINGOLD (D WI)'], 'FEINGOLD (D WI)', col = rgb(0, 0, 1, .5), cex = .8)
	points(articles$unity[articles$senator == 'MCCAIN (R AZ)'], articles$mentions[articles$senator == 'MCCAIN (R AZ)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'SNOWE (R ME)'], articles$mentions[articles$senator == 'SNOWE (R ME)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'STEVENS (R AK)'], articles$mentions[articles$senator == 'STEVENS (R AK)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'MIKULSKI (D MD)'], articles$mentions[articles$senator == 'MIKULSKI (D MD)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'LIEBERMAN (D CT)'], articles$mentions[articles$senator == 'LIEBERMAN (D CT)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'BIDEN (D DE)'], articles$mentions[articles$senator == 'BIDEN (D DE)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'SANTORUM (R PA)'], articles$mentions[articles$senator == 'SANTORUM (R PA)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'SPECTER (R PA)'], articles$mentions[articles$senator == 'SPECTER (R PA)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'DOMENICI (R NM)'], articles$mentions[articles$senator == 'DOMENICI (R NM)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	points(articles$unity[articles$senator == 'FEINGOLD (D WI)'], articles$mentions[articles$senator == 'FEINGOLD (D WI)'], pch = 16, col = rgb(1, 0, 0, .55), cex = .8)
	abline(lm(articles$mentions ~ articles$unity), col = rgb(0, 1, 0, .25), lwd = 3)
dev.off()

## main graphic replication ##
	data <- read.dta('dataCompleteWithBetas.dta')
	data1 <- data
	data <- data[c(1:5000), c(1:87)]
		
## grab mean values ##

	data1 <- data1[is.na(data1$trueYea) == FALSE &
		is.na(data1$trueNay) == FALSE &
		is.na(data1$partyYea) == FALSE &
		is.na(data1$polInt) == FALSE &
		is.na(data1$preferYea) == FALSE &
		is.na(data1$preferNay) == FALSE &
		is.na(data1$partyAgree) == FALSE &
		is.na(data1$female) == FALSE &
		is.na(data1$race) == FALSE &
		is.na(data1$income) == FALSE &
		is.na(data1$edu) == FALSE &
		is.na(data1$abort) == FALSE &
		is.na(data1$gains) == FALSE &
		is.na(data1$imm) == FALSE &
		is.na(data1$iraq) == FALSE &
		is.na(data1$stem) == FALSE &
		is.na(data1$wage) == FALSE, ]

## grab mean covariate values for graphics ##

	yes <- mean(data1$yes, na.rm = TRUE)
	no <- mean(data1$no, na.rm = TRUE)
	trueYea <- mean(data1$trueYea, na.rm = TRUE)
	trueNay <- mean(data1$trueNay, na.rm = TRUE)
	partyYea <- mean(data1$partyYea, na.rm = TRUE)
	pastUnity <- mean(data1$pastUnity, na.rm = TRUE)
	freshman <- mean(data1$freshman, na.rm = TRUE)
	polInt <- mean(data1$polInt, na.rm = TRUE)
	preferYea <- mean(data1$preferYea, na.rm = TRUE)
	preferNay <- mean(data1$preferNay, na.rm = TRUE)
	partyAgree <- mean(data1$partyAgree, na.rm = TRUE)
	female <- mean(data1$female, na.rm = TRUE)
	race <- mean(data1$race, na.rm = TRUE)
	income <- mean(data1$income, na.rm = TRUE)
	edu <- mean(data1$edu, na.rm = TRUE)
	abort <- mean(data1$abort, na.rm = TRUE)
	gains <- mean(data1$gains, na.rm = TRUE)
	imm <- mean(data1$imm, na.rm = TRUE)
	iraq <- mean(data1$iraq, na.rm = TRUE)
	stem <- mean(data1$stem, na.rm = TRUE)
	wage <- mean(data1$wage, na.rm = TRUE)
	
## 2x2, true x party effects over party unity ##

	range <- seq(0.75, 1.00, 0.005)
	
	yesA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))


## full model pics ##

	for(i in 1:length(range)){

	yes <- data$bbt1 * 0 +
		data$bbt3 * 1 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 0 * range[i] + 
		data$bbt15 * 0 * 0 + 
		data$bbt17 * 1 * range[i] + 
		data$bbt19 * 1 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 0 * polInt + 
		data$bbt31 * 1 * polInt + 
		data$bbt33 * 1 * polInt + 
		data$bbt35 * 0 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 1 * range[i] * polInt + 
		data$bbt41 * 1 * 0 * polInt + 
		data$bbt43 * 1 * range[i] * polInt + 
		data$bbt45 * 1 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 0 * 0.5 + 
		data$bbt55 * 1 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
	no <- data$bbt2 * 0 + 
		data$bbt4 * 1 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 0 * range[i] + 
		data$bbt16 * 0 * 0 + 
		data$bbt18 * 1 * range[i] + 
		data$bbt20 * 1 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 0 * polInt + 
		data$bbt32 * 1 * polInt + 
		data$bbt34 * 1 * polInt + 
		data$bbt36 * 0 * range[i] * polInt + 
		data$bbt38 * 0 * 0 * polInt + 
		data$bbt40 * 1 * range[i] * polInt + 
		data$bbt42 * 1 * 0 * polInt + 
		data$bbt44 * 1 * range[i] * polInt + 
		data$bbt46 * 1 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 0 * 0.5 + 
		data$bbt56 * 1 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87
		
	yesA[, i] <- yes
	noA[, i] <- no
	
	denom <- 1 + exp(yesA[, i]) + exp(noA[, i])
	yesA[, i] <- exp(yesA[, i]) / denom
	noA[, i] <- exp(noA[, i]) / denom
	dkA[, i] <- 1 - yesA[, i] - noA[, i]

	##########***********##########
	
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 1 * polInt + 
		data$bbt31 * 0 * polInt + 
		data$bbt33 * 1 * polInt + 
		data$bbt35 * 1 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 0 * range[i] * polInt + 
		data$bbt41 * 0 * 0 * polInt + 
		data$bbt43 * 1 * range[i] * polInt + 
		data$bbt45 * 1 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 1 * polInt + 
		data$bbt32 * 0 * polInt + 
		data$bbt34 * 1 * polInt + 
		data$bbt36 * 1 * range[i] * polInt + 
		data$bbt38 * 1 * 0 * polInt + 
		data$bbt40 * 0 * range[i] * polInt + 
		data$bbt42 * 0 * 0 * polInt + 
		data$bbt44 * 1 * range[i] * polInt + 
		data$bbt46 * 1 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87
		
	yesB[, i] <- yes
	noB[, i] <- no
	
	denom <- 1 + exp(yesB[, i]) + exp(noB[, i])
	yesB[, i] <- exp(yesB[, i]) / denom
	noB[, i] <- exp(noB[, i]) / denom
	dkB[, i] <- 1 - yesB[, i] - noB[, i]

	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 1 * polInt + 
		data$bbt31 * 0 * polInt + 
		data$bbt33 * 0 * polInt + 
		data$bbt35 * 1 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 0 * range[i] * polInt + 
		data$bbt41 * 0 * 0 * polInt + 
		data$bbt43 * 0 * range[i] * polInt + 
		data$bbt45 * 0 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 1 * polInt + 
		data$bbt32 * 0 * polInt + 
		data$bbt34 * 0 * polInt + 
		data$bbt36 * 1 * range[i] * polInt + 
		data$bbt38 * 1 * 0 * polInt + 
		data$bbt40 * 0 * range[i] * polInt + 
		data$bbt42 * 0 * 0 * polInt + 
		data$bbt44 * 0 * range[i] * polInt + 
		data$bbt46 * 0 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87

		
	yesC[, i] <- yes
	noC[, i] <- no
	
	denom <- 1 + exp(yesC[, i]) + exp(noC[, i])
	yesC[, i] <- exp(yesC[, i]) / denom
	noC[, i] <- exp(noC[, i]) / denom
	dkC[, i] <- 1 - yesC[, i] - noC[, i]

	##########***********##########
	
	
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 1 * polInt + 
		data$bbt31 * 0 * polInt + 
		data$bbt33 * 1 * polInt + 
		data$bbt35 * 1 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 0 * range[i] * polInt + 
		data$bbt41 * 0 * 0 * polInt + 
		data$bbt43 * 1 * range[i] * polInt + 
		data$bbt45 * 1 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 1 * polInt + 
		data$bbt32 * 0 * polInt + 
		data$bbt34 * 1 * polInt + 
		data$bbt36 * 1 * range[i] * polInt + 
		data$bbt38 * 1 * 0 * polInt + 
		data$bbt40 * 0 * range[i] * polInt + 
		data$bbt42 * 0 * 0 * polInt + 
		data$bbt44 * 1 * range[i] * polInt + 
		data$bbt46 * 1 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87

	yesD[, i] <- yes
	noD[, i] <- no
	
	denom <- 1 + exp(yesD[, i]) + exp(noD[, i])
	yesD[, i] <- exp(yesD[, i]) / denom
	noD[, i] <- exp(noD[, i]) / denom
	dkD[, i] <- 1 - yesD[, i] - noD[, i]
		
	}

pdf('figure2.pdf', width = 10, height = 12)
par(mfrow = c(2, 1))

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = 'Holding True Vote at Yea and\nChanging Party Line from Nay to Yea',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesD[, i]-yesC[, i], 0.975)
		hiN[i] <- quantile(noD[, i]-noC[, i], 0.975)
		hiD[i] <- quantile(dkD[, i]-dkC[, i], 0.975)
		loY[i] <- quantile(yesD[, i]-yesC[, i], 0.025)
		loN[i] <- quantile(noD[, i]-noC[, i], 0.025)
		loD[i] <- quantile(dkD[, i]-dkC[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesD-yesC), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noD-noC), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkD-dkC), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkD-dkC)[51], 'DK ')
	text(1.01, colMeans(noD-noC)[51], 'Nay')
	text(1.01, colMeans(yesD-yesC)[51], 'Yea')
	text(0.74, colMeans(dkD-dkC)[1], 'DK ')
	text(0.74, colMeans(noD-noC)[1], 'Nay')
	text(0.74, colMeans(yesD-yesC)[1], 'Yea')

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = 'Holding Party Line at Yea and\nChanging True Vote from Nay to Yea',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesB[, i]-yesA[, i], 0.975)
		hiN[i] <- quantile(noB[, i]-noA[, i], 0.975)
		hiD[i] <- quantile(dkB[, i]-dkA[, i], 0.975)
		loY[i] <- quantile(yesB[, i]-yesA[, i], 0.025)
		loN[i] <- quantile(noB[, i]-noA[, i], 0.025)
		loD[i] <- quantile(dkB[, i]-dkA[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesB-yesA), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noB-noA), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkB-dkA), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkB-dkA)[51], 'DK ')
	text(1.01, colMeans(noB-noA)[51], 'Nay')
	text(1.01, colMeans(yesB-yesA)[51], 'Yea')
	text(0.74, colMeans(dkB-dkA)[1], 'DK ')
	text(0.74, colMeans(noB-noA)[1], 'Nay')
	text(0.74, colMeans(yesB-yesA)[1], 'Yea')

dev.off()

## now figure 3 ##

	for(i in 1:length(range)){

	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 1 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * 1 + 
		data$bbt27 * 0 * 1 + 
		data$bbt29 * 1 * 1 + 
		data$bbt31 * 0 * 1 + 
		data$bbt33 * 0 * 1 + 
		data$bbt35 * 1 * range[i] * 1 + 
		data$bbt37 * 1 * 0 * 1 + 
		data$bbt39 * 0 * range[i] * 1 + 
		data$bbt41 * 0 * 0 * 1 + 
		data$bbt43 * 0 * range[i] * 1 + 
		data$bbt45 * 0 * 0 * 1 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 1 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * 1 + 
		data$bbt28 * 0 * 1 + 
		data$bbt30 * 1 * 1 + 
		data$bbt32 * 0 * 1 + 
		data$bbt34 * 0 * 1 + 
		data$bbt36 * 1 * range[i] * 1 + 
		data$bbt38 * 1 * 0 * 1 + 
		data$bbt40 * 0 * range[i] * 1 + 
		data$bbt42 * 0 * 0 * 1 + 
		data$bbt44 * 0 * range[i] * 1 + 
		data$bbt46 * 0 * 0 * 1 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87

		
	yesA[, i] <- yes
	noA[, i] <- no
	
	denom <- 1 + exp(yesA[, i]) + exp(noA[, i])
	yesA[, i] <- exp(yesA[, i]) / denom
	noA[, i] <- exp(noA[, i]) / denom
	dkA[, i] <- 1 - yesA[, i] - noA[, i]

	##########***********##########
	
	
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 1 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * 1 + 
		data$bbt27 * 0 * 1 + 
		data$bbt29 * 1 * 1 + 
		data$bbt31 * 0 * 1 + 
		data$bbt33 * 1 * 1 + 
		data$bbt35 * 1 * range[i] * 1 + 
		data$bbt37 * 1 * 0 * 1 + 
		data$bbt39 * 0 * range[i] * 1 + 
		data$bbt41 * 0 * 0 * 1 + 
		data$bbt43 * 1 * range[i] * 1 + 
		data$bbt45 * 1 * 0 * 1 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 1 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * 1 + 
		data$bbt28 * 0 * 1 + 
		data$bbt30 * 1 * 1 + 
		data$bbt32 * 0 * 1 + 
		data$bbt34 * 1 * 1 + 
		data$bbt36 * 1 * range[i] * 1 + 
		data$bbt38 * 1 * 0 * 1 + 
		data$bbt40 * 0 * range[i] * 1 + 
		data$bbt42 * 0 * 0 * 1 + 
		data$bbt44 * 1 * range[i] * 1 + 
		data$bbt46 * 1 * 0 * 1 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87

	yesB[, i] <- yes
	noB[, i] <- no
	
	denom <- 1 + exp(yesB[, i]) + exp(noB[, i])
	yesB[, i] <- exp(yesB[, i]) / denom
	noB[, i] <- exp(noB[, i]) / denom
	dkB[, i] <- 1 - yesB[, i] - noB[, i]


	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 3 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * 3 + 
		data$bbt27 * 0 * 3 + 
		data$bbt29 * 1 * 3 + 
		data$bbt31 * 0 * 3 + 
		data$bbt33 * 0 * 3 + 
		data$bbt35 * 1 * range[i] * 3 + 
		data$bbt37 * 1 * 0 * 3 + 
		data$bbt39 * 0 * range[i] * 3 + 
		data$bbt41 * 0 * 0 * 3 + 
		data$bbt43 * 0 * range[i] * 3 + 
		data$bbt45 * 0 * 0 * 3 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 3 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * 3 + 
		data$bbt28 * 0 * 3 + 
		data$bbt30 * 1 * 3 + 
		data$bbt32 * 0 * 3 + 
		data$bbt34 * 0 * 3 + 
		data$bbt36 * 1 * range[i] * 3 + 
		data$bbt38 * 1 * 0 * 3 + 
		data$bbt40 * 0 * range[i] * 3 + 
		data$bbt42 * 0 * 0 * 3 + 
		data$bbt44 * 0 * range[i] * 3 + 
		data$bbt46 * 0 * 0 * 3 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87

		
	yesC[, i] <- yes
	noC[, i] <- no
	
	denom <- 1 + exp(yesC[, i]) + exp(noC[, i])
	yesC[, i] <- exp(yesC[, i]) / denom
	noC[, i] <- exp(noC[, i]) / denom
	dkC[, i] <- 1 - yesC[, i] - noC[, i]

	##########***********##########
	
	
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 3 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * 3 + 
		data$bbt27 * 0 * 3 + 
		data$bbt29 * 1 * 3 + 
		data$bbt31 * 0 * 3 + 
		data$bbt33 * 1 * 3 + 
		data$bbt35 * 1 * range[i] * 3 + 
		data$bbt37 * 1 * 0 * 3 + 
		data$bbt39 * 0 * range[i] * 3 + 
		data$bbt41 * 0 * 0 * 3 + 
		data$bbt43 * 1 * range[i] * 3 + 
		data$bbt45 * 1 * 0 * 3 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 3 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * 3 + 
		data$bbt28 * 0 * 3 + 
		data$bbt30 * 1 * 3 + 
		data$bbt32 * 0 * 3 + 
		data$bbt34 * 1 * 3 + 
		data$bbt36 * 1 * range[i] * 3 + 
		data$bbt38 * 1 * 0 * 3 + 
		data$bbt40 * 0 * range[i] * 3 + 
		data$bbt42 * 0 * 0 * 3 + 
		data$bbt44 * 1 * range[i] * 3 + 
		data$bbt46 * 1 * 0 * 3 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87

	yesD[, i] <- yes
	noD[, i] <- no
	
	denom <- 1 + exp(yesD[, i]) + exp(noD[, i])
	yesD[, i] <- exp(yesD[, i]) / denom
	noD[, i] <- exp(noD[, i]) / denom
	dkD[, i] <- 1 - yesD[, i] - noD[, i]
		
	} 


pdf('figure3.pdf', width = 10, height = 12)
par(mfrow = c(2, 1))

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = '',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesB[, i]-yesA[, i], 0.975)
		hiN[i] <- quantile(noB[, i]-noA[, i], 0.975)
		hiD[i] <- quantile(dkB[, i]-dkA[, i], 0.975)
		loY[i] <- quantile(yesB[, i]-yesA[, i], 0.025)
		loN[i] <- quantile(noB[, i]-noA[, i], 0.025)
		loD[i] <- quantile(dkB[, i]-dkA[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesB-yesA), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noB-noA), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkB-dkA), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkB-dkA)[51], 'DK ')
	text(1.01, colMeans(noB-noA)[51], 'Nay')
	text(1.01, colMeans(yesB-yesA)[51], 'Yea')
	text(0.74, colMeans(dkB-dkA)[1]+0.025, 'DK ')
	text(0.74, colMeans(noB-noA)[1]-0.025, 'Nay')
	text(0.74, colMeans(yesB-yesA)[1], 'Yea')
	mtext('Low Interest Voters', line = 2.5, cex = 1.5, font = 2)
	mtext('Holding True Vote at Yea and\nChanging Party Line from Nay to Yea', line = 0.25)

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = '',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesD[, i]-yesC[, i], 0.975)
		hiN[i] <- quantile(noD[, i]-noC[, i], 0.975)
		hiD[i] <- quantile(dkD[, i]-dkC[, i], 0.975)
		loY[i] <- quantile(yesD[, i]-yesC[, i], 0.025)
		loN[i] <- quantile(noD[, i]-noC[, i], 0.025)
		loD[i] <- quantile(dkD[, i]-dkC[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesD-yesC), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noD-noC), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkD-dkC), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkD-dkC)[51], 'DK ')
	text(1.01, colMeans(noD-noC)[51], 'Nay')
	text(1.01, colMeans(yesD-yesC)[51], 'Yea')
	text(0.74, colMeans(dkD-dkC)[1], 'DK ')
	text(0.74, colMeans(noD-noC)[1], 'Nay')
	text(0.74, colMeans(yesD-yesC)[1], 'Yea')
	mtext('High Interest Voters', line = 2.5, cex = 1.5, font = 2)
	mtext('Holding True Vote at Yea and\nChanging Party Line from Nay to Yea', line = 0.25)
dev.off()

## figure 4 ##
## Interest Effects ##
## this will change interest from low to high ##
## for the case of a party line matching the true vote ##
## as well as for a mistmatch ##

	for(i in 1:length(range)){

	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 1 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * 1 + 
		data$bbt27 * 0 * 1 + 
		data$bbt29 * 1 * 1 + 
		data$bbt31 * 0 * 1 + 
		data$bbt33 * 1 * 1 + 
		data$bbt35 * 1 * range[i] * 1 + 
		data$bbt37 * 1 * 0 * 1 + 
		data$bbt39 * 0 * range[i] * 1 + 
		data$bbt41 * 0 * 0 * 1 + 
		data$bbt43 * 1 * range[i] * 1 + 
		data$bbt45 * 1 * 0 * 1 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 1 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * 1 + 
		data$bbt28 * 0 * 1 + 
		data$bbt30 * 1 * 1 + 
		data$bbt32 * 0 * 1 + 
		data$bbt34 * 1 * 1 + 
		data$bbt36 * 1 * range[i] * 1 + 
		data$bbt38 * 1 * 0 * 1 + 
		data$bbt40 * 0 * range[i] * 1 + 
		data$bbt42 * 0 * 0 * 1 + 
		data$bbt44 * 1 * range[i] * 1 + 
		data$bbt46 * 1 * 0 * 1 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87
		
	yesA[, i] <- yes
	noA[, i] <- no
	
	denom <- 1 + exp(yesA[, i]) + exp(noA[, i])
	yesA[, i] <- exp(yesA[, i]) / denom
	noA[, i] <- exp(noA[, i]) / denom
	dkA[, i] <- 1 - yesA[, i] - noA[, i]

	##########***********##########
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 3 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * 3 + 
		data$bbt27 * 0 * 3 + 
		data$bbt29 * 1 * 3 + 
		data$bbt31 * 0 * 3 + 
		data$bbt33 * 1 * 3 + 
		data$bbt35 * 1 * range[i] * 3 + 
		data$bbt37 * 1 * 0 * 3 + 
		data$bbt39 * 0 * range[i] * 3 + 
		data$bbt41 * 0 * 0 * 3 + 
		data$bbt43 * 1 * range[i] * 3 + 
		data$bbt45 * 1 * 0 * 3 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 3 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * 3 + 
		data$bbt28 * 0 * 3 + 
		data$bbt30 * 1 * 3 + 
		data$bbt32 * 0 * 3 + 
		data$bbt34 * 1 * 3 + 
		data$bbt36 * 1 * range[i] * 3 + 
		data$bbt38 * 1 * 0 * 3 + 
		data$bbt40 * 0 * range[i] * 3 + 
		data$bbt42 * 0 * 0 * 3 + 
		data$bbt44 * 1 * range[i] * 3 + 
		data$bbt46 * 1 * 0 * 3 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87
		
	yesB[, i] <- yes
	noB[, i] <- no
	
	denom <- 1 + exp(yesB[, i]) + exp(noB[, i])
	yesB[, i] <- exp(yesB[, i]) / denom
	noB[, i] <- exp(noB[, i]) / denom
	dkB[, i] <- 1 - yesB[, i] - noB[, i]
	
		
	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 1 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * 1 + 
		data$bbt27 * 0 * 1 + 
		data$bbt29 * 1 * 1 + 
		data$bbt31 * 0 * 1 + 
		data$bbt33 * 0 * 1 + 
		data$bbt35 * 1 * range[i] * 1 + 
		data$bbt37 * 1 * 0 * 1 + 
		data$bbt39 * 0 * range[i] * 1 + 
		data$bbt41 * 0 * 0 * 1 + 
		data$bbt43 * 0 * range[i] * 1 + 
		data$bbt45 * 0 * 0 * 1 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 1 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * 1 + 
		data$bbt28 * 0 * 1 + 
		data$bbt30 * 1 * 1 + 
		data$bbt32 * 0 * 1 + 
		data$bbt34 * 0 * 1 + 
		data$bbt36 * 1 * range[i] * 1 + 
		data$bbt38 * 1 * 0 * 1 + 
		data$bbt40 * 0 * range[i] * 1 + 
		data$bbt42 * 0 * 0 * 1 + 
		data$bbt44 * 0 * range[i] * 1 + 
		data$bbt46 * 0 * 0 * 1 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87
	
	yesC[, i] <- yes
	noC[, i] <- no
	
	denom <- 1 + exp(yesC[, i]) + exp(noC[, i])
	yesC[, i] <- exp(yesC[, i]) / denom
	noC[, i] <- exp(noC[, i]) / denom
	dkC[, i] <- 1 - yesC[, i] - noC[, i]
		
		
	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 3 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * 3 + 
		data$bbt27 * 0 * 3 + 
		data$bbt29 * 1 * 3 + 
		data$bbt31 * 0 * 3 + 
		data$bbt33 * 0 * 3 + 
		data$bbt35 * 1 * range[i] * 3 + 
		data$bbt37 * 1 * 0 * 3 + 
		data$bbt39 * 0 * range[i] * 3 + 
		data$bbt41 * 0 * 0 * 3 + 
		data$bbt43 * 0 * range[i] * 3 + 
		data$bbt45 * 0 * 0 * 3 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83 + 
		data$bbt85
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 3 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * 3 + 
		data$bbt28 * 0 * 3 + 
		data$bbt30 * 1 * 3 + 
		data$bbt32 * 0 * 3 + 
		data$bbt34 * 0 * 3 + 
		data$bbt36 * 1 * range[i] * 3 + 
		data$bbt38 * 1 * 0 * 3 + 
		data$bbt40 * 0 * range[i] * 3 + 
		data$bbt42 * 0 * 0 * 3 + 
		data$bbt44 * 0 * range[i] * 3 + 
		data$bbt46 * 0 * 0 * 3 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84 + 
		data$bbt86 + 
		data$bbt87
	
	yesD[, i] <- yes
	noD[, i] <- no
	
	denom <- 1 + exp(yesD[, i]) + exp(noD[, i])
	yesD[, i] <- exp(yesD[, i]) / denom
	noD[, i] <- exp(noD[, i]) / denom
	dkD[, i] <- 1 - yesD[, i] - noD[, i]
		
	}
	
## create the plots ##

pdf('figure4.pdf', width = 10, height = 12)
par(mfrow = c(2, 1))

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = '',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesB[, i]-yesA[, i], 0.975)
		hiN[i] <- quantile(noB[, i]-noA[, i], 0.975)
		hiD[i] <- quantile(dkB[, i]-dkA[, i], 0.975)
		loY[i] <- quantile(yesB[, i]-yesA[, i], 0.025)
		loN[i] <- quantile(noB[, i]-noA[, i], 0.025)
		loD[i] <- quantile(dkB[, i]-dkA[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesB-yesA), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noB-noA), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkB-dkA), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkB-dkA)[51], 'DK ')
	text(1.01, colMeans(noB-noA)[51], 'Nay')
	text(1.01, colMeans(yesB-yesA)[51], 'Yea')
	text(0.74, colMeans(dkB-dkA)[1]+0.025, 'DK ')
	text(0.74, colMeans(noB-noA)[1]-0.025, 'Nay')
	text(0.74, colMeans(yesB-yesA)[1], 'Yea')
	mtext('Effect of Political Interest', line = 2.5, cex = 1.5, font = 2)
	mtext('Holding True Vote at Yea and Party Line at Yea\nChanging Political Interest from Low to High', line = 0.25)

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = '',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesD[, i]-yesC[, i], 0.975)
		hiN[i] <- quantile(noD[, i]-noC[, i], 0.975)
		hiD[i] <- quantile(dkD[, i]-dkC[, i], 0.975)
		loY[i] <- quantile(yesD[, i]-yesC[, i], 0.025)
		loN[i] <- quantile(noD[, i]-noC[, i], 0.025)
		loD[i] <- quantile(dkD[, i]-dkC[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesD-yesC), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noD-noC), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkD-dkC), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkD-dkC)[51], 'DK ')
	text(1.01, colMeans(noD-noC)[51], 'Nay')
	text(1.01, colMeans(yesD-yesC)[51], 'Yea')
	text(0.74, colMeans(dkD-dkC)[1], 'DK ')
	text(0.74, colMeans(noD-noC)[1], 'Nay')
	text(0.74, colMeans(yesD-yesC)[1], 'Yea')
	mtext('Effect of Political Interest', line = 2.5, cex = 1.5, font = 2)
	mtext('Holding True Vote at Yea and Party Line at Nay\nChanging Political Interest from Low to High', line = 0.25)
dev.off()


## triplots ##
## figure 5 ##
	
	A <- c()
	B <- c()
	C <- c()
	
	for(i in 1:5000){
	
		data2 <- data1[rank(runif(length(data1$polInt), 1, length(data1$polInt)))[c(1:30)], ]
				
		yes <- data$bbt1[i] * 0.5 + 
			data$bbt3[i] * 0.5 + 
			data$bbt5[i] * 1 + 
			data$bbt7[i] * 0.75 + 
			data$bbt9[i] * 0 + 
			data$bbt11[i] * data2$polInt + 
			data$bbt13[i] * 0.5 * 0.75 + 
			data$bbt15[i] * 0.5 * 0 + 
			data$bbt17[i] * 0.5 * 0.75 + 
			data$bbt19[i] * 0.5 * 0 + 
			data$bbt21[i] * 1 * 0.75 + 
			data$bbt23[i] * 1 * 0 + 
			data$bbt25[i] * 0.75 * data2$polInt + 
			data$bbt27[i] * 0 * data2$polInt + 
			data$bbt29[i] * 0.5 * data2$polInt + 
			data$bbt31[i] * 0.5 * data2$polInt + 
			data$bbt33[i] * 1 * data2$polInt + 
			data$bbt35[i] * 0.5 * 0.75 * data2$polInt + 
			data$bbt37[i] * 0.5 * 0 * data2$polInt + 
			data$bbt39[i] * 0.5 * 0.75 * data2$polInt + 
			data$bbt41[i] * 0.5 * 0 * data2$polInt + 
			data$bbt43[i] * 1 * 0.75 * data2$polInt + 
			data$bbt45[i] * 1 * 0 * data2$polInt + 
			data$bbt47[i] * data2$preferYea + 
			data$bbt49[i] * data2$preferNay + 
			data$bbt51[i] * data2$partyAgree + 
			data$bbt53[i] * 0.5 * data2$partyAgree + 
			data$bbt55[i] * 0.5 * data2$partyAgree + 
			data$bbt57[i] * data2$partyYea * data2$partyAgree + 
			data$bbt59[i] * data2$preferYea * data2$partyAgree + 
			data$bbt61[i] * data2$preferNay * data2$partyAgree + 
			data$bbt63[i] * data2$female + 
			data$bbt65[i] * data2$race + 
			data$bbt67[i] * data2$income + 
			data$bbt69[i] * data2$edu + 
			data$bbt71[i] * data2$abort + 
			data$bbt73[i] * data2$gains + 
			data$bbt75[i] * data2$imm + 
			data$bbt77[i] * data2$iraq + 
			data$bbt79[i] * data2$stem + 
			data$bbt81[i] * data2$wage + 
			data$bbt83[i] +
			data$bbt85[i]
	
	
		no <- data$bbt2[i] * 0.5 + 
			data$bbt4[i] * 0.5 + 
			data$bbt6[i] * 1 + 
			data$bbt8[i] * 0.75 + 
			data$bbt10[i] * 0 + 
			data$bbt12[i] * data2$polInt + 
			data$bbt14[i] * 0.5 * 0.75 + 
			data$bbt16[i] * 0.5 * 0 + 
			data$bbt18[i] * 0.5 * 0.75 + 
			data$bbt20[i] * 0.5 * 0 + 
			data$bbt22[i] * 1 * 0.75 + 
			data$bbt24[i] * 1 * 0 + 
			data$bbt26[i] * 0.75 * data2$polInt + 
			data$bbt28[i] * 0 * data2$polInt + 
			data$bbt30[i] * 0.5 * data2$polInt + 
			data$bbt32[i] * 0.5 * data2$polInt + 
			data$bbt34[i] * 1 * data2$polInt + 
			data$bbt36[i] * 0.5 * 0.75 * data2$polInt + 
			data$bbt38[i] * 0.5 * 0 * data2$polInt + 
			data$bbt40[i] * 0.5 * 0.75 * data2$polInt + 
			data$bbt42[i] * 0.5 * 0 * data2$polInt + 
			data$bbt44[i] * 1 * 0.75 * data2$polInt + 
			data$bbt46[i] * 1 * 0 * data2$polInt + 
			data$bbt48[i] * data2$preferYea + 
			data$bbt50[i] * data2$preferNay + 
			data$bbt52[i] * data2$partyAgree + 
			data$bbt54[i] * 0.5 * data2$partyAgree + 
			data$bbt56[i] * 0.5 * data2$partyAgree + 
			data$bbt58[i] * data2$partyYea * data2$partyAgree + 
			data$bbt60[i] * data2$preferYea * data2$partyAgree + 
			data$bbt62[i] * data2$preferNay * data2$partyAgree + 
			data$bbt64[i] * data2$female + 
			data$bbt66[i] * data2$race + 
			data$bbt68[i] * data2$income + 
			data$bbt70[i] * data2$edu + 
			data$bbt72[i] * data2$abort + 
			data$bbt74[i] * data2$gains + 
			data$bbt76[i] * data2$imm + 
			data$bbt78[i] * data2$iraq + 
			data$bbt80[i] * data2$stem + 
			data$bbt82[i] * data2$wage + 
			data$bbt84[i] +
			data$bbt86[i] + 
			data$bbt87[i]
		
	
		yesA <- yes
		noA <- no
		
		denom <- 1 + exp(yesA) + exp(noA)
		yesA <- exp(yesA) / denom
		noA <- exp(noA) / denom
		dkA <- 1 - yesA - noA
		
		a <- cbind(noA, dkA, yesA)
		A <- rbind(A, a)
		
		yes <- data$bbt1[i] * 0.5 + 
			data$bbt3[i] * 0.5 + 
			data$bbt5[i] * 1 + 
			data$bbt7[i] * 0.95 + 
			data$bbt9[i] * 0 + 
			data$bbt11[i] * data2$polInt + 
			data$bbt13[i] * 0.5 * 0.95 + 
			data$bbt15[i] * 0.5 * 0 + 
			data$bbt17[i] * 0.5 * 0.95 + 
			data$bbt19[i] * 0.5 * 0 + 
			data$bbt21[i] * 1 * 0.95 + 
			data$bbt23[i] * 1 * 0 + 
			data$bbt25[i] * 0.95 * data2$polInt + 
			data$bbt27[i] * 0 * data2$polInt + 
			data$bbt29[i] * 0.5 * data2$polInt + 
			data$bbt31[i] * 0.5 * data2$polInt + 
			data$bbt33[i] * 1 * data2$polInt + 
			data$bbt35[i] * 0.5 * 0.95 * data2$polInt + 
			data$bbt37[i] * 0.5 * 0 * data2$polInt + 
			data$bbt39[i] * 0.5 * 0.95 * data2$polInt + 
			data$bbt41[i] * 0.5 * 0 * data2$polInt + 
			data$bbt43[i] * 1 * 0.95 * data2$polInt + 
			data$bbt45[i] * 1 * 0 * data2$polInt + 
			data$bbt47[i] * data2$preferYea + 
			data$bbt49[i] * data2$preferNay + 
			data$bbt51[i] * data2$partyAgree + 
			data$bbt53[i] * 0.5 * data2$partyAgree + 
			data$bbt55[i] * 0.5 * data2$partyAgree + 
			data$bbt57[i] * data2$partyYea * data2$partyAgree + 
			data$bbt59[i] * data2$preferYea * data2$partyAgree + 
			data$bbt61[i] * data2$preferNay * data2$partyAgree + 
			data$bbt63[i] * data2$female + 
			data$bbt65[i] * data2$race + 
			data$bbt67[i] * data2$income + 
			data$bbt69[i] * data2$edu + 
			data$bbt71[i] * data2$abort + 
			data$bbt73[i] * data2$gains + 
			data$bbt75[i] * data2$imm + 
			data$bbt77[i] * data2$iraq + 
			data$bbt79[i] * data2$stem + 
			data$bbt81[i] * data2$wage + 
			data$bbt83[i] +
			data$bbt85[i]
	
	
		no <- data$bbt2[i] * 0.5 + 
			data$bbt4[i] * 0.5 + 
			data$bbt6[i] * 1 + 
			data$bbt8[i] * 0.95 + 
			data$bbt10[i] * 0 + 
			data$bbt12[i] * data2$polInt + 
			data$bbt14[i] * 0.5 * 0.95 + 
			data$bbt16[i] * 0.5 * 0 + 
			data$bbt18[i] * 0.5 * 0.95 + 
			data$bbt20[i] * 0.5 * 0 + 
			data$bbt22[i] * 1 * 0.95 + 
			data$bbt24[i] * 1 * 0 + 
			data$bbt26[i] * 0.95 * data2$polInt + 
			data$bbt28[i] * 0 * data2$polInt + 
			data$bbt30[i] * 0.5 * data2$polInt + 
			data$bbt32[i] * 0.5 * data2$polInt + 
			data$bbt34[i] * 1 * data2$polInt + 
			data$bbt36[i] * 0.5 * 0.95 * data2$polInt + 
			data$bbt38[i] * 0.5 * 0 * data2$polInt + 
			data$bbt40[i] * 0.5 * 0.95 * data2$polInt + 
			data$bbt42[i] * 0.5 * 0 * data2$polInt + 
			data$bbt44[i] * 1 * 0.95 * data2$polInt + 
			data$bbt46[i] * 1 * 0 * data2$polInt + 
			data$bbt48[i] * data2$preferYea + 
			data$bbt50[i] * data2$preferNay + 
			data$bbt52[i] * data2$partyAgree + 
			data$bbt54[i] * 0.5 * data2$partyAgree + 
			data$bbt56[i] * 0.5 * data2$partyAgree + 
			data$bbt58[i] * data2$partyYea * data2$partyAgree + 
			data$bbt60[i] * data2$preferYea * data2$partyAgree + 
			data$bbt62[i] * data2$preferNay * data2$partyAgree + 
			data$bbt64[i] * data2$female + 
			data$bbt66[i] * data2$race + 
			data$bbt68[i] * data2$income + 
			data$bbt70[i] * data2$edu + 
			data$bbt72[i] * data2$abort + 
			data$bbt74[i] * data2$gains + 
			data$bbt76[i] * data2$imm + 
			data$bbt78[i] * data2$iraq + 
			data$bbt80[i] * data2$stem + 
			data$bbt82[i] * data2$wage + 
			data$bbt84[i] +
			data$bbt86[i] + 
			data$bbt87[i]
		
	
		yesB <- yes
		noB <- no
		
		denom <- 1 + exp(yesB) + exp(noB)
		yesB <- exp(yesB) / denom
		noB <- exp(noB) / denom
		dkB <- 1 - yesB - noB
		
		b <- cbind(noB, dkB, yesB)
		B <- rbind(B, b)
		
		yes <- data$bbt1[i] * 0.5 + 
			data$bbt3[i] * 0.5 + 
			data$bbt5[i] * 1 + 
			data$bbt7[i] * 0 + 
			data$bbt9[i] * 1 + 
			data$bbt11[i] * data2$polInt + 
			data$bbt13[i] * 0.5 * 0 + 
			data$bbt15[i] * 0.5 * 1 + 
			data$bbt17[i] * 0.5 * 0 + 
			data$bbt19[i] * 0.5 * 1 + 
			data$bbt21[i] * 1 * 0 + 
			data$bbt23[i] * 1 * 1 + 
			data$bbt25[i] * 0 * data2$polInt + 
			data$bbt27[i] * 1 * data2$polInt + 
			data$bbt29[i] * 0.5 * data2$polInt + 
			data$bbt31[i] * 0.5 * data2$polInt + 
			data$bbt33[i] * 1 * data2$polInt + 
			data$bbt35[i] * 0.5 * 0 * data2$polInt + 
			data$bbt37[i] * 0.5 * 1 * data2$polInt + 
			data$bbt39[i] * 0.5 * 0 * data2$polInt + 
			data$bbt41[i] * 0.5 * 1 * data2$polInt + 
			data$bbt43[i] * 1 * 0 * data2$polInt + 
			data$bbt45[i] * 1 * 1 * data2$polInt + 
			data$bbt47[i] * data2$preferYea + 
			data$bbt49[i] * data2$preferNay + 
			data$bbt51[i] * data2$partyAgree + 
			data$bbt53[i] * 0.5 * data2$partyAgree + 
			data$bbt55[i] * 0.5 * data2$partyAgree + 
			data$bbt57[i] * data2$partyYea * data2$partyAgree + 
			data$bbt59[i] * data2$preferYea * data2$partyAgree + 
			data$bbt61[i] * data2$preferNay * data2$partyAgree + 
			data$bbt63[i] * data2$female + 
			data$bbt65[i] * data2$race + 
			data$bbt67[i] * data2$income + 
			data$bbt69[i] * data2$edu + 
			data$bbt71[i] * data2$abort + 
			data$bbt73[i] * data2$gains + 
			data$bbt75[i] * data2$imm + 
			data$bbt77[i] * data2$iraq + 
			data$bbt79[i] * data2$stem + 
			data$bbt81[i] * data2$wage + 
			data$bbt83[i] +
			data$bbt85[i]
			
	
		no <- data$bbt2[i] * 0.5 + 
			data$bbt4[i] * 0.5 + 
			data$bbt6[i] * 1 + 
			data$bbt8[i] * 0 + 
			data$bbt10[i] * 1 + 
			data$bbt12[i] * data2$polInt + 
			data$bbt14[i] * 0.5 * 0 + 
			data$bbt16[i] * 0.5 * 1 + 
			data$bbt18[i] * 0.5 * 0 + 
			data$bbt20[i] * 0.5 * 1 + 
			data$bbt22[i] * 1 * 0 + 
			data$bbt24[i] * 1 * 1 + 
			data$bbt26[i] * 0 * data2$polInt + 
			data$bbt28[i] * 1 * data2$polInt + 
			data$bbt30[i] * 0.5 * data2$polInt + 
			data$bbt32[i] * 0.5 * data2$polInt + 
			data$bbt34[i] * 1 * data2$polInt + 
			data$bbt36[i] * 0.5 * 0 * data2$polInt + 
			data$bbt38[i] * 0.5 * 1 * data2$polInt + 
			data$bbt40[i] * 0.5 * 0 * data2$polInt + 
			data$bbt42[i] * 0.5 * 1 * data2$polInt + 
			data$bbt44[i] * 1 * 0 * data2$polInt + 
			data$bbt46[i] * 1 * 1 * data2$polInt + 
			data$bbt48[i] * data2$preferYea + 
			data$bbt50[i] * data2$preferNay + 
			data$bbt52[i] * data2$partyAgree + 
			data$bbt54[i] * 0.5 * data2$partyAgree + 
			data$bbt56[i] * 0.5 * data2$partyAgree + 
			data$bbt58[i] * data2$partyYea * data2$partyAgree + 
			data$bbt60[i] * data2$preferYea * data2$partyAgree + 
			data$bbt62[i] * data2$preferNay * data2$partyAgree + 
			data$bbt64[i] * data2$female + 
			data$bbt66[i] * data2$race + 
			data$bbt68[i] * data2$income + 
			data$bbt70[i] * data2$edu + 
			data$bbt72[i] * data2$abort + 
			data$bbt74[i] * data2$gains + 
			data$bbt76[i] * data2$imm + 
			data$bbt78[i] * data2$iraq + 
			data$bbt80[i] * data2$stem + 
			data$bbt82[i] * data2$wage + 
			data$bbt84[i] +
			data$bbt86[i] + 
			data$bbt87[i]
			
	
		yesC <- yes
		noC <- no
		
		denom <- 1 + exp(yesC) + exp(noC)
		yesC <- exp(yesC) / denom
		noC <- exp(noC) / denom
		dkC <- 1 - yesC - noC
		
		c <- cbind(noC, dkC, yesC)
		C <- rbind(C, c)
	
	
		cat(i, 'complete...\n')	
	}

	A <- A[rank(runif(length(A[, 1]), 1, length(A[, 1])))[1:20000], ]
	B <- B[rank(runif(length(B[, 1]), 1, length(B[, 1])))[1:20000], ]
	C <- C[rank(runif(length(C[, 1]), 1, length(C[, 1])))[1:20000], ]

## the triplot() command can be a little wonky, so we'll just hard code it here ##

triplot <- function(x, y=NULL, z=NULL, 
         labels=dimnames(x)[[2]], 
         txt=dimnames(x)[[1]], legend=NULL,
         legend.split=NULL,
         inner=TRUE, inner.col=c('lightblue','pink'),
         inner.lty=c(2,3), add=FALSE, main="", ...){

  old.par <- par(xpd=TRUE)
  on.exit(par(old.par))

  if( is.data.frame(x) ) x <- as.matrix(x)
  x <- cbind(x,y,z)
  if( ncol(x) < 2 || ncol(x) > 3 ){
    stop("need 2 or 3 columns")
  }
  if( ncol(x)==3 ){
    x <- sweep(x,1,FUN="/",apply(x,1,sum))
  }
  if( ncol(x)==2 ){
    x <- cbind(x, 1-x[,1]-x[,2])
  }

  if(dev.cur()==1){
    win.graph()
    add <- FALSE
  }
	
  if( !add ){
    
    xstar <-  -0.03327397
    
    plot( c(0,1,2,0), c(0,sqrt(3),0,0), type="l",
         lwd=1.5, xlim=c(-xstar,2+xstar),
         xlab="",ylab="",axes=FALSE, main=main)
	
    if(inner){
      lines( c(1,1.5,0.5,1), c(0,sqrt(3)/2,sqrt(3)/2,0),
            lwd=.5, col=inner.col[1], lty=inner.lty[1])
      lines( c(1.25, 1, .75, 1.25), 
            c(sqrt(3)/4, sqrt(3)/2, sqrt(3)/4, sqrt(3)/4),
            lwd=0.25, col=inner.col[2],lty=inner.lty[2])
    }
	
    if(length(labels)==0){
      labels <- c("X","Y","Z")
    }
	
    ystar <- 0.0797549
    text( c(0,2,1), c(-ystar,-ystar,sqrt(3)+(ystar/2)),
         labels, cex=1.25) 
  }
	
  newy <- x[,3] * sqrt(3)
  newx <- 2-2*x[,1]-x[,3]
	
  if(length(txt)==length(newx)){
    text(newx,newy,txt,...)
  } else {
    points(newx,newy,...)
  }
	
  if(length(legend)==length(newx)){
    labpos <- function(y){
      strh <- par("cxy")[2]*1.15
      y2 <- sort(y)
      df <- y2[-1] - y2[-length(y2)]
      i <- 1
      while(any (df < strh)){
        y2[c(df < strh, FALSE)] <- y2[ c(df < strh,FALSE)] - strh/10
        y2[c(FALSE, df < strh)] <- y2[ c(FALSE,df < strh)] + strh/10
        if(min(y2)<0){y2 <- y2 - min(y2)}
        y2 <- sort(y2)
        df <- y2[-1] - y2[ -length(y2)]
        i <- i+1
        if(i>100){break}
      }
      y2
    }

    if(length(legend.split)==1){
      tmp.x <- quantile(newx, legend.split)
      y1 <- newy[newx <= tmp.x]
      y1 <- labpos(y1)[order(order(y1))]
      text(rep(-0.01,length(y1)), y1, legend[newx<=tmp.x], adj=1)
      segments(rep(0,length(y1)), y1, newx[newx<=tmp.x], newy[newx<=tmp.x])

      y2 <- newy[newx>tmp.x]
      y2 <- labpos(y2)[order(order(y2))]
      text(rep(2.01,length(y2)), y2, legend[newx>tmp.x], adj=0)
      segments(rep(2,length(y2)), y2, newx[newx>tmp.x], newy[newx>tmp.x])
    } else {
      if(any(newx <= 1)){
        y1 <- newy[newx<=1]
        y1 <- labpos(y1)[order(order(y1))]
        text(rep(-0.01,length(y1)), y1, legend[newx<=1],adj=1)
        segments(rep(0,length(y1)), y1, newx[newx<=1], newy[newx<=1])
      }
      if(any(newx > 1)){
        y2 <- newy[newx>1]
        y2 <- labpos(y2)[order(order(y2))]
        text(rep(2.01,length(y2)), y2, legend[newx>1],adj=0)
        segments(rep(2,length(y2)), y2, newx[newx>1], newy[newx>1])
      }
    }
  }
  
  invisible(cbind(x=newx,y=newy))
}

pdf('figure5.pdf', height = 5, width = 15)
	par(mfrow = c(1, 3))
	triplot(A[, c(1:3)],
		col = rgb(1, 0, 0, 0.005), cex = 2, pch = 16,
		labels = c('Against Party', 'DK', 'With Party'), inner.col = rgb(0, 0, 0, 0.5))
		mtext('Mavericks', cex = 1.1, font = 2, line = 1)


	triplot(B[, c(1:3)],
		col = rgb(1, 0, 0, 0.005), cex = 2, pch = 16,
		labels = c('Against Party', 'DK', 'With Party'), inner.col = rgb(0, 0, 0, 0.5))
		mtext('Loyal Partisans', cex = 1.1, font = 2, line = 1)

	triplot(C[, c(1:3)],
		col = rgb(1, 0, 0, 0.005), cex = 2, pch = 16,
		labels = c('Against Party', 'DK', 'With Party'), inner.col = rgb(0, 0, 0, 0.5))
		mtext('Freshmen', cex = 1.1, font = 2, line = 1)
dev.off()


##########################################################################################
##################  this portion produces the graphics in the appendix  ##################
##########################################################################################

## appendix figure a1.1 ##

	data <- read.dta('dataComplete.dta')
	state <- c()
	sen1 <- c()
	sen2 <- c()
	unity1 <- c()
	unity2 <- c()
	ratio1 <- c()
	ratio2 <- c()
	sig <- c()
	
	for(i in c(1:50)){
		state[i] <- unique(data$state)[i]
		sen1[i] <- unique(data$senator[data$state == unique(data$state)[i]])[1]
		sen2[i] <- unique(data$senator[data$state == unique(data$state)[i]])[2]
		unity1[i] <- unique(data$pastUnity[data$senator == sen1[i]])
		unity2[i] <- unique(data$pastUnity[data$senator == sen2[i]])
		ratio1[i] <- length(data$apply[data$senator == sen1[i] & data$apply == 1])/length(data$apply[data$senator == sen1[i] & data$apply == 0])
		ratio2[i] <- length(data$apply[data$senator == sen2[i] & data$apply == 1])/length(data$apply[data$senator == sen2[i] & data$apply == 0])
	}
	
	state <- state[-12]
	sen1 <- sen1[-12]
	sen2 <- sen2[-12]
	unity1 <- unity1[-12]
	unity2 <- unity2[-12]
	ratio1 <- ratio1[-12]
	ratio2 <- ratio2[-12]

	pdf('perceivedLoyaltyComparison.pdf')
	plot(unity1-unity2, ratio1-ratio2,
		xlim = c(-.25, .25),
		ylim = c(-1.5, 1.5),
		xlab = 'Difference in Unity Scores',
		ylab = 'Difference in Perceived Loyalty',
		main = 'Comparing the Relationship Between Party Unity\n and Perceived Loyalty Within States')
	lines(c(-.5, .5), c(0, 0))
	lines(c(0, 0), c(-1.5, 1.5))
	dev.off()
	
## appendix figures a3.1 and a3.2 ##
	
	ses <- read.table('seClusterResults.csv', skip = 4, sep = ',')
	pull <- seq(2, 144, by = 2)
	names <- ses[-pull, 1]
	ses <- ses[pull, ]

	pdf('a3.1.pdf', width = 9, height = 6)
	par(mfrow = c(2, 3))
	hist(((ses$V3-ses$V2)/ses$V2)*100,
		main = 'SEs Clustered on: Senator',
		xlab = 'Percent Change in SE',
		col = rgb(1, 0, 0, 0.5),
		border = rgb(1, 0, 0, 0.25))

	hist(((ses$V4-ses$V2)/ses$V2)*100,
		main = 'SEs Clustered on: Issue',
		xlab = 'Percent Change in SE',
		col = rgb(1, 0, 0, 0.5),
		border = rgb(1, 0, 0, 0.25))

	hist(((ses$V5-ses$V2)/ses$V2)*100,
		main = 'SEs Clustered on: Respondent',
		xlab = 'Percent Change in SE',
		col = rgb(1, 0, 0, 0.5),
		border = rgb(1, 0, 0, 0.25))

	hist(((ses$V6-ses$V2)/ses$V2)*100,
		main = 'SEs Clustered on: State',
		xlab = 'Percent Change in SE',
		col = rgb(1, 0, 0, 0.5),
		border = rgb(1, 0, 0, 0.25))

	hist(((ses$V7-ses$V2)/ses$V2)*100,
		main = 'SEs Clustered on: Senator-Issue',
		xlab = 'Percent Change in SE',
		col = rgb(1, 0, 0, 0.5),
		border = rgb(1, 0, 0, 0.25))	
	dev.off()
	
## now a3.2	##
	
	no <- ses[seq(2, 72, 2), ]
	yes <- ses[seq(1, 71, 2), ]
	
	pdf('a3.2.pdf', width = 8, height = 8)
		par(mfrow = c(2, 1))
		plot(1, 1, type = 'n',
			xlim = c(1, 36),
			ylim = c(-200, 200),
			axes = FALSE,
			ylab = 'Percent Change in SE',
			xlab = '')
		par(las = 2, cex.axis = 0.5)
		axis(2, at = seq(-100, 200, 50))
		axis(1, at = 1:36, labels = names[seq(1, 71, 2)], line = -3)
		points(1:36, ((yes$V3-yes$V2)/yes$V2)*100, col = rgb(1, 0, 0, 0.5), pch = 16)
		lines(1:36, ((yes$V3-yes$V2)/yes$V2)*100, col = rgb(1, 0, 0, 0.5))

		points(1:36, ((yes$V5-yes$V2)/yes$V2)*100, col = rgb(0, 1, 0, 0.5), pch = 16)
		lines(1:36, ((yes$V5-yes$V2)/yes$V2)*100, col = rgb(0, 1, 0, 0.5))

		points(1:36, ((yes$V6-yes$V2)/yes$V2)*100, col = rgb(0, 0, 1, 0.5), pch = 16)
		lines(1:36, ((yes$V6-yes$V2)/yes$V2)*100, col = rgb(0, 0, 1, 0.5))

		points(1:36, ((yes$V7-yes$V2)/yes$V2)*100, col = rgb(1, 1, 0, 0.75), pch = 16)
		lines(1:36, ((yes$V7-yes$V2)/yes$V2)*100, col = rgb(1, 1, 0, 0.75))
		legend('topleft', legend = c('Senator', 'Respondent', 'State', 'Senator-Issue'),
			fill = c(rgb(1, 0, 0, 0.5), rgb(0, 1, 0, 0.5), rgb(0, 0, 1, 0.5), rgb(1, 1, 0, 0.75)),
			border = FALSE, bty = 'n')

		plot(1, 1, type = 'n',
			xlim = c(1, 36),
			ylim = c(-200, 200),
			axes = FALSE,
			ylab = 'Percent Change in SE',
			xlab = '')
		axis(2, at = seq(-100, 200, 50))
		axis(1, at = 1:36, labels = names[seq(2, 72, 2)], line = -3)
		points(1:36, ((no$V3-no$V2)/no$V2)*100, col = rgb(1, 0, 0, 0.5), pch = 16)
		lines(1:36, ((no$V3-no$V2)/no$V2)*100, col = rgb(1, 0, 0, 0.5))

		points(1:36, ((no$V5-no$V2)/no$V2)*100, col = rgb(0, 1, 0, 0.5), pch = 16)
		lines(1:36, ((no$V5-no$V2)/no$V2)*100, col = rgb(0, 1, 0, 0.5))

		points(1:36, ((no$V6-no$V2)/no$V2)*100, col = rgb(0, 0, 1, 0.5), pch = 16)
		lines(1:36, ((no$V6-no$V2)/no$V2)*100, col = rgb(0, 0, 1, 0.5))

		points(1:36, ((no$V7-no$V2)/no$V2)*100, col = rgb(1, 1, 0, 0.75), pch = 16)
		lines(1:36, ((no$V7-no$V2)/no$V2)*100, col = rgb(1, 1, 0, 0.75))
		legend('topleft', legend = c('Senator', 'Respondent', 'State', 'Senator-Issue'),
			fill = c(rgb(1, 0, 0, 0.5), rgb(0, 1, 0, 0.5), rgb(0, 0, 1, 0.5), rgb(1, 1, 0, 0.75)),
			border = FALSE, bty = 'n')
		dev.off()
			
## appendix figures a4.1 ##

	data <- read.dta('dataCompleteWithBetas.dta')
	dataX <- read.dta('simpleModelBetas.dta')
	dataX <- dataX[c(1:5000), c(1:28)]

	data1 <- data
	data <- data[c(1:5000), c(1:87)]
		
## grab mean values ##

	data1 <- data1[is.na(data1$trueYea) == FALSE &
		is.na(data1$trueNay) == FALSE &
		is.na(data1$partyYea) == FALSE &
		is.na(data1$polInt) == FALSE &
		is.na(data1$preferYea) == FALSE &
		is.na(data1$preferNay) == FALSE &
		is.na(data1$partyAgree) == FALSE &
		is.na(data1$female) == FALSE &
		is.na(data1$race) == FALSE &
		is.na(data1$income) == FALSE &
		is.na(data1$edu) == FALSE &
		is.na(data1$abort) == FALSE &
		is.na(data1$gains) == FALSE &
		is.na(data1$imm) == FALSE &
		is.na(data1$iraq) == FALSE &
		is.na(data1$stem) == FALSE &
		is.na(data1$wage) == FALSE, ]

	names(data)
	yes <- mean(data1$yes, na.rm = TRUE)
	no <- mean(data1$no, na.rm = TRUE)
	trueYea <- mean(data1$trueYea, na.rm = TRUE)
	trueNay <- mean(data1$trueNay, na.rm = TRUE)
	partyYea <- mean(data1$partyYea, na.rm = TRUE)
	pastUnity <- mean(data1$pastUnity, na.rm = TRUE)
	freshman <- mean(data1$freshman, na.rm = TRUE)
	polInt <- mean(data1$polInt, na.rm = TRUE)
	preferYea <- mean(data1$preferYea, na.rm = TRUE)
	preferNay <- mean(data1$preferNay, na.rm = TRUE)
	partyAgree <- mean(data1$partyAgree, na.rm = TRUE)
	female <- mean(data1$female, na.rm = TRUE)
	race <- mean(data1$race, na.rm = TRUE)
	income <- mean(data1$income, na.rm = TRUE)
	edu <- mean(data1$edu, na.rm = TRUE)
	abort <- mean(data1$abort, na.rm = TRUE)
	gains <- mean(data1$gains, na.rm = TRUE)
	imm <- mean(data1$imm, na.rm = TRUE)
	iraq <- mean(data1$iraq, na.rm = TRUE)
	stem <- mean(data1$stem, na.rm = TRUE)
	wage <- mean(data1$wage, na.rm = TRUE)
	
## 2x2, true x party effects over party unity ##

	range <- seq(0.75, 1.00, 0.005)
	
	yesA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

## probability change pics ##
## simple model ##
	for(i in 1:length(range)){

	yes <- dataX$bbt1 * 0 +
		dataX$bbt3 * 1 + 
		dataX$bbt5 * 1 + 
		dataX$bbt7 * range[i] + 
		dataX$bbt9 * 0 + 
		dataX$bbt11 * 0 * range[i] + 
		dataX$bbt13 * 0 * 0 + 
		dataX$bbt15 * 1 * range[i] + 
		dataX$bbt17 * 1 * 0 + 
		dataX$bbt19 * 1 * range[i] + 
		dataX$bbt21 * 1 * 0 + 
		dataX$bbt23 * 0.5 + 
		dataX$bbt25 * 0.5 + 
		dataX$bbt27
		
	no <- dataX$bbt2 * 0 +
		dataX$bbt4 * 1 + 
		dataX$bbt6 * 1 + 
		dataX$bbt8 * range[i] + 
		dataX$bbt10 * 0 + 
		dataX$bbt12 * 0 * range[i] + 
		dataX$bbt14 * 0 * 0 + 
		dataX$bbt16 * 1 * range[i] + 
		dataX$bbt18 * 1 * 0 + 
		dataX$bbt20 * 1 * range[i] + 
		dataX$bbt22 * 1 * 0 + 
		dataX$bbt24 * 0.5 + 
		dataX$bbt26 * 0.5 + 
		dataX$bbt28
		
	yesA[, i] <- yes
	noA[, i] <- no
	
	denom <- 1 + exp(yesA[, i]) + exp(noA[, i])
	yesA[, i] <- exp(yesA[, i]) / denom
	noA[, i] <- exp(noA[, i]) / denom
	dkA[, i] <- 1 - yesA[, i] - noA[, i]

	##########***********##########
	
	yes <- dataX$bbt1 * 1 +
		dataX$bbt3 * 0 + 
		dataX$bbt5 * 1 + 
		dataX$bbt7 * range[i] + 
		dataX$bbt9 * 0 + 
		dataX$bbt11 * 1 * range[i] + 
		dataX$bbt13 * 1 * 0 + 
		dataX$bbt15 * 0 * range[i] + 
		dataX$bbt17 * 0 * 0 + 
		dataX$bbt19 * 1 * range[i] + 
		dataX$bbt21 * 1 * 0 + 
		dataX$bbt23 * 0.5 + 
		dataX$bbt25 * 0.5 + 
		dataX$bbt27
		
	no <- dataX$bbt2 * 1 +
		dataX$bbt4 * 0 + 
		dataX$bbt6 * 1 + 
		dataX$bbt8 * range[i] + 
		dataX$bbt10 * 0 + 
		dataX$bbt12 * 1 * range[i] + 
		dataX$bbt14 * 1 * 0 + 
		dataX$bbt16 * 0 * range[i] + 
		dataX$bbt18 * 0 * 0 + 
		dataX$bbt20 * 1 * range[i] + 
		dataX$bbt22 * 1 * 0 + 
		dataX$bbt24 * 0.5 + 
		dataX$bbt26 * 0.5 + 
		dataX$bbt28
		
	yesB[, i] <- yes
	noB[, i] <- no
	
	denom <- 1 + exp(yesB[, i]) + exp(noB[, i])
	yesB[, i] <- exp(yesB[, i]) / denom
	noB[, i] <- exp(noB[, i]) / denom
	dkB[, i] <- 1 - yesB[, i] - noB[, i]

	yes <- dataX$bbt1 * 1 +
		dataX$bbt3 * 0 + 
		dataX$bbt5 * 0 + 
		dataX$bbt7 * range[i] + 
		dataX$bbt9 * 0 + 
		dataX$bbt11 * 1 * range[i] + 
		dataX$bbt13 * 1 * 0 + 
		dataX$bbt15 * 0 * range[i] + 
		dataX$bbt17 * 0 * 0 + 
		dataX$bbt19 * 0 * range[i] + 
		dataX$bbt21 * 0 * 0 + 
		dataX$bbt23 * 0.5 + 
		dataX$bbt25 * 0.5 + 
		dataX$bbt27
		
	no <- dataX$bbt2 * 1 +
		dataX$bbt4 * 0 + 
		dataX$bbt6 * 0 + 
		dataX$bbt8 * range[i] + 
		dataX$bbt10 * 0 + 
		dataX$bbt12 * 1 * range[i] + 
		dataX$bbt14 * 1 * 0 + 
		dataX$bbt16 * 0 * range[i] + 
		dataX$bbt18 * 0 * 0 + 
		dataX$bbt20 * 0 * range[i] + 
		dataX$bbt22 * 0 * 0 + 
		dataX$bbt24 * 0.5 + 
		dataX$bbt26 * 0.5 + 
		dataX$bbt28

		
	yesC[, i] <- yes
	noC[, i] <- no
	
	denom <- 1 + exp(yesC[, i]) + exp(noC[, i])
	yesC[, i] <- exp(yesC[, i]) / denom
	noC[, i] <- exp(noC[, i]) / denom
	dkC[, i] <- 1 - yesC[, i] - noC[, i]

	##########***********##########
	
	
	
	yes <- dataX$bbt1 * 1 +
		dataX$bbt3 * 0 + 
		dataX$bbt5 * 1 + 
		dataX$bbt7 * range[i] + 
		dataX$bbt9 * 0 + 
		dataX$bbt11 * 1 * range[i] + 
		dataX$bbt13 * 1 * 0 + 
		dataX$bbt15 * 0 * range[i] + 
		dataX$bbt17 * 0 * 0 + 
		dataX$bbt19 * 1 * range[i] + 
		dataX$bbt21 * 1 * 0 + 
		dataX$bbt23 * 0.5 + 
		dataX$bbt25 * 0.5 + 
		dataX$bbt27
		
	no <- dataX$bbt2 * 1 +
		dataX$bbt4 * 0 + 
		dataX$bbt6 * 1 + 
		dataX$bbt8 * range[i] + 
		dataX$bbt10 * 0 + 
		dataX$bbt12 * 1 * range[i] + 
		dataX$bbt14 * 1 * 0 + 
		dataX$bbt16 * 0 * range[i] + 
		dataX$bbt18 * 0 * 0 + 
		dataX$bbt20 * 1 * range[i] + 
		dataX$bbt22 * 1 * 0 + 
		dataX$bbt24 * 0.5 + 
		dataX$bbt26 * 0.5 + 
		dataX$bbt28

	yesD[, i] <- yes
	noD[, i] <- no
	
	denom <- 1 + exp(yesD[, i]) + exp(noD[, i])
	yesD[, i] <- exp(yesD[, i]) / denom
	noD[, i] <- exp(noD[, i]) / denom
	dkD[, i] <- 1 - yesD[, i] - noD[, i]
		
	}

pdf('a4.1.pdf', width = 10, height = 12)
par(mfrow = c(2, 1))

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = 'Holding True Vote at Yea and\nChanging Party Line from Nay to Yea',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesD[, i]-yesC[, i], 0.975)
		hiN[i] <- quantile(noD[, i]-noC[, i], 0.975)
		hiD[i] <- quantile(dkD[, i]-dkC[, i], 0.975)
		loY[i] <- quantile(yesD[, i]-yesC[, i], 0.025)
		loN[i] <- quantile(noD[, i]-noC[, i], 0.025)
		loD[i] <- quantile(dkD[, i]-dkC[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesD-yesC), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noD-noC), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkD-dkC), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkD-dkC)[51], 'DK ')
	text(1.01, colMeans(noD-noC)[51], 'Nay')
	text(1.01, colMeans(yesD-yesC)[51], 'Yea')
	text(0.74, colMeans(dkD-dkC)[1], 'DK ')
	text(0.74, colMeans(noD-noC)[1], 'Nay')
	text(0.74, colMeans(yesD-yesC)[1]-0.02, 'Yea')

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Response Probability',
		main = 'Holding Party Line at Yea and\nChanging True Vote from Nay to Yea',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesB[, i]-yesA[, i], 0.975)
		hiN[i] <- quantile(noB[, i]-noA[, i], 0.975)
		hiD[i] <- quantile(dkB[, i]-dkA[, i], 0.975)
		loY[i] <- quantile(yesB[, i]-yesA[, i], 0.025)
		loN[i] <- quantile(noB[, i]-noA[, i], 0.025)
		loD[i] <- quantile(dkB[, i]-dkA[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesB-yesA), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noB-noA), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkB-dkA), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkB-dkA)[51]-0.01, 'DK ')
	text(1.01, colMeans(noB-noA)[51]+0.01, 'Nay')
	text(1.01, colMeans(yesB-yesA)[51], 'Yea')
	text(0.74, colMeans(dkB-dkA)[1], 'DK ')
	text(0.74, colMeans(noB-noA)[1], 'Nay')
	text(0.74, colMeans(yesB-yesA)[1], 'Yea')
dev.off()



## appendix figures a4.2 ##

	data <- read.dta('dataMaverickCompleteWithBetas.dta')
	summary(data)
	data1 <- data
	data <- data[c(1:5000), c(1:71)]
		
## grab mean values ##

	data1 <- data1[is.na(data1$trueYea) == FALSE &
		is.na(data1$trueNay) == FALSE &
		is.na(data1$partyYea) == FALSE &
		is.na(data1$polInt) == FALSE &
		is.na(data1$preferYea) == FALSE &
		is.na(data1$preferNay) == FALSE &
		is.na(data1$maverick) == FALSE &
		is.na(data1$partyAgree) == FALSE &
		is.na(data1$female) == FALSE &
		is.na(data1$race) == FALSE &
		is.na(data1$income) == FALSE &
		is.na(data1$edu) == FALSE &
		is.na(data1$abort) == FALSE &
		is.na(data1$gains) == FALSE &
		is.na(data1$imm) == FALSE &
		is.na(data1$iraq) == FALSE &
		is.na(data1$stem) == FALSE &
		is.na(data1$wage) == FALSE &
		data1$freshman == 0, ]

	names(data)
	yes <- mean(data1$yes, na.rm = TRUE)
	no <- mean(data1$no, na.rm = TRUE)
	trueYea <- mean(data1$trueYea, na.rm = TRUE)
	trueNay <- mean(data1$trueNay, na.rm = TRUE)
	partyYea <- mean(data1$partyYea, na.rm = TRUE)
	maverick <- mean(data1$maverick, na.rm = TRUE)
	polInt <- mean(data1$polInt, na.rm = TRUE)
	preferYea <- mean(data1$preferYea, na.rm = TRUE)
	preferNay <- mean(data1$preferNay, na.rm = TRUE)
	partyAgree <- mean(data1$partyAgree, na.rm = TRUE)
	female <- mean(data1$female, na.rm = TRUE)
	race <- mean(data1$race, na.rm = TRUE)
	income <- mean(data1$income, na.rm = TRUE)
	edu <- mean(data1$edu, na.rm = TRUE)
	abort <- mean(data1$abort, na.rm = TRUE)
	gains <- mean(data1$gains, na.rm = TRUE)
	imm <- mean(data1$imm, na.rm = TRUE)
	iraq <- mean(data1$iraq, na.rm = TRUE)
	stem <- mean(data1$stem, na.rm = TRUE)
	wage <- mean(data1$wage, na.rm = TRUE)
	
## 2x2, true x party effects over party unity ##

	range <- seq(0, 0.13, 0.001)
	
	yesA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	for(i in 1:length(range)){

	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * polInt + 
		data$bbt11 * 1 * range[i] + 
		data$bbt13 * 0 * range[i] + 
		data$bbt15 * 0 * range[i] + 
		data$bbt17 * range[i] * polInt + 
		data$bbt19 * 1 * polInt + 
		data$bbt21 * 0 * polInt + 
		data$bbt23 * 0 * polInt + 
		data$bbt25 * 1 * range[i] * polInt + 
		data$bbt27 * 0 * range[i] * polInt + 
		data$bbt29 * 0 * range[i] * polInt + 
		data$bbt31 * 0.5 + 
		data$bbt33 * 0.5 + 
		data$bbt35 * 0.5 + 
		data$bbt37 * 1 * 0.5 + 
		data$bbt39 * 0 * 0.5 + 
		data$bbt41 * 0 * 0.5 + 
		data$bbt43 * 0.5 * 0.5 + 
		data$bbt45 * 0.5 * 0.5 + 
		data$bbt47 * female + 
		data$bbt49 * race + 
		data$bbt51 * income + 
		data$bbt53 * edu + 
		data$bbt55 * abort + 
		data$bbt57 * gains + 
		data$bbt59 * imm + 
		data$bbt61 * iraq + 
		data$bbt63 * stem + 
		data$bbt65 * wage + 
		data$bbt67 + 
		data$bbt69
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * polInt + 
		data$bbt12 * 1 * range[i] + 
		data$bbt14 * 0 * range[i] + 
		data$bbt16 * 0 * range[i] + 
		data$bbt18 * range[i] * polInt + 
		data$bbt20 * 1 * polInt + 
		data$bbt22 * 0 * polInt + 
		data$bbt24 * 0 * polInt + 
		data$bbt26 * 1 * range[i] * polInt + 
		data$bbt28 * 0 * range[i] * polInt + 
		data$bbt30 * 0 * range[i] * polInt + 
		data$bbt32 * 0.5 + 
		data$bbt34 * 0.5 + 
		data$bbt36 * 0.5 + 
		data$bbt38 * 1 * 0.5 + 
		data$bbt40 * 0 * 0.5 + 
		data$bbt42 * 0 * 0.5 + 
		data$bbt44 * 0.5 * 0.5 + 
		data$bbt46 * 0.5 * 0.5 + 
		data$bbt48 * female + 
		data$bbt50 * race + 
		data$bbt52 * income + 
		data$bbt54 * edu + 
		data$bbt56 * abort + 
		data$bbt58 * gains + 
		data$bbt60 * imm + 
		data$bbt62 * iraq + 
		data$bbt64 * stem + 
		data$bbt66 * wage + 
		data$bbt68 + 
		data$bbt70 +
		data$bbt71
		
	yesA[, i] <- yes
	noA[, i] <- no
	
	denom <- 1 + exp(yesA[, i]) + exp(noA[, i])
	yesA[, i] <- exp(yesA[, i]) / denom
	noA[, i] <- exp(noA[, i]) / denom
	dkA[, i] <- 1 - yesA[, i] - noA[, i]

	##########***********##########
	
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * polInt + 
		data$bbt11 * 1 * range[i] + 
		data$bbt13 * 0 * range[i] + 
		data$bbt15 * 1 * range[i] + 
		data$bbt17 * range[i] * polInt + 
		data$bbt19 * 1 * polInt + 
		data$bbt21 * 0 * polInt + 
		data$bbt23 * 1 * polInt + 
		data$bbt25 * 1 * range[i] * polInt + 
		data$bbt27 * 0 * range[i] * polInt + 
		data$bbt29 * 1 * range[i] * polInt + 
		data$bbt31 * 0.5 + 
		data$bbt33 * 0.5 + 
		data$bbt35 * 0.5 + 
		data$bbt37 * 1 * 0.5 + 
		data$bbt39 * 0 * 0.5 + 
		data$bbt41 * 1 * 0.5 + 
		data$bbt43 * 0.5 * 0.5 + 
		data$bbt45 * 0.5 * 0.5 + 
		data$bbt47 * female + 
		data$bbt49 * race + 
		data$bbt51 * income + 
		data$bbt53 * edu + 
		data$bbt55 * abort + 
		data$bbt57 * gains + 
		data$bbt59 * imm + 
		data$bbt61 * iraq + 
		data$bbt63 * stem + 
		data$bbt65 * wage + 
		data$bbt67 + 
		data$bbt69
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * polInt + 
		data$bbt12 * 1 * range[i] + 
		data$bbt14 * 0 * range[i] + 
		data$bbt16 * 1 * range[i] + 
		data$bbt18 * range[i] * polInt + 
		data$bbt20 * 1 * polInt + 
		data$bbt22 * 0 * polInt + 
		data$bbt24 * 1 * polInt + 
		data$bbt26 * 1 * range[i] * polInt + 
		data$bbt28 * 0 * range[i] * polInt + 
		data$bbt30 * 1 * range[i] * polInt + 
		data$bbt32 * 0.5 + 
		data$bbt34 * 0.5 + 
		data$bbt36 * 0.5 + 
		data$bbt38 * 1 * 0.5 + 
		data$bbt40 * 0 * 0.5 + 
		data$bbt42 * 1 * 0.5 + 
		data$bbt44 * 0.5 * 0.5 + 
		data$bbt46 * 0.5 * 0.5 + 
		data$bbt48 * female + 
		data$bbt50 * race + 
		data$bbt52 * income + 
		data$bbt54 * edu + 
		data$bbt56 * abort + 
		data$bbt58 * gains + 
		data$bbt60 * imm + 
		data$bbt62 * iraq + 
		data$bbt64 * stem + 
		data$bbt66 * wage + 
		data$bbt68 + 
		data$bbt70 +
		data$bbt71		
				
	yesB[, i] <- yes
	noB[, i] <- no
	
	denom <- 1 + exp(yesB[, i]) + exp(noB[, i])
	yesB[, i] <- exp(yesB[, i]) / denom
	noB[, i] <- exp(noB[, i]) / denom
	dkB[, i] <- 1 - yesB[, i] - noB[, i]
		
	}


pdf('a4.2.pdf', width = 10, height = 6)
	plot(1, 1, type = 'n',
		xlab = 'Proportion of Media Messages Regarding Maverick Behavior',
		ylab = 'Change in Predicted Response Probability',
		main = 'Holding True Vote at Yea and\nChanging Party Line from Nay to Yea',
		xlim = c(-0.005, 0.135),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0, 0.13, 0.01))
	axis(2)
	polygon(c(0, 0, 0.13, 0.13), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0, 0.13), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 0:13){
		lines(c(0+(0.01*i), 0+(0.01*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesB[, i]-yesA[, i], 0.975)
		hiN[i] <- quantile(noB[, i]-noA[, i], 0.975)
		hiD[i] <- quantile(dkB[, i]-dkA[, i], 0.975)
		loY[i] <- quantile(yesB[, i]-yesA[, i], 0.025)
		loN[i] <- quantile(noB[, i]-noA[, i], 0.025)
		loD[i] <- quantile(dkB[, i]-dkA[, i], 0.025)
	}
	lines(c(0, 0.13), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesB-yesA), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noB-noA), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkB-dkA), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(0.135, colMeans(dkB-dkA)[131], 'DK ')
	text(0.135, colMeans(noB-noA)[131], 'Nay')
	text(0.135, colMeans(yesB-yesA)[131], 'Yea')
	text(-0.005, colMeans(dkB-dkA)[1], 'DK ')
	text(-0.005, colMeans(noB-noA)[1], 'Nay')
	text(-0.005, colMeans(yesB-yesA)[1], 'Yea')
dev.off()

## now the full sample figures a4.3 and a4.4 ##

	data <- read.dta('dataCompleteWithBetasCorrIssueDumAllFULLSAMPLE.dta')

	data1 <- data
	data <- data[c(1:5000), c(1:87)]
		
## grab mean values ##

	data1 <- data1[is.na(data1$trueYea) == FALSE &
		is.na(data1$trueNay) == FALSE &
		is.na(data1$partyYea) == FALSE &
		is.na(data1$polInt) == FALSE &
		is.na(data1$preferYea) == FALSE &
		is.na(data1$preferNay) == FALSE &
		is.na(data1$partyAgree) == FALSE &
		is.na(data1$female) == FALSE &
		is.na(data1$race) == FALSE &
		is.na(data1$income) == FALSE &
		is.na(data1$edu) == FALSE &
		is.na(data1$abort) == FALSE &
		is.na(data1$gains) == FALSE &
		is.na(data1$imm) == FALSE &
		is.na(data1$iraq) == FALSE &
		is.na(data1$stem) == FALSE &
		is.na(data1$wage) == FALSE, ]

	names(data)
	yes <- mean(data1$yes, na.rm = TRUE)
	no <- mean(data1$no, na.rm = TRUE)
	trueYea <- mean(data1$trueYea, na.rm = TRUE)
	trueNay <- mean(data1$trueNay, na.rm = TRUE)
	partyYea <- mean(data1$partyYea, na.rm = TRUE)
	pastUnity <- mean(data1$pastUnity, na.rm = TRUE)
	freshman <- mean(data1$freshman, na.rm = TRUE)
	polInt <- mean(data1$polInt, na.rm = TRUE)
	preferYea <- mean(data1$preferYea, na.rm = TRUE)
	preferNay <- mean(data1$preferNay, na.rm = TRUE)
	partyAgree <- mean(data1$partyAgree, na.rm = TRUE)
	female <- mean(data1$female, na.rm = TRUE)
	race <- mean(data1$race, na.rm = TRUE)
	income <- mean(data1$income, na.rm = TRUE)
	edu <- mean(data1$edu, na.rm = TRUE)
	abort <- mean(data1$abort, na.rm = TRUE)
	gains <- mean(data1$gains, na.rm = TRUE)
	imm <- mean(data1$imm, na.rm = TRUE)
	iraq <- mean(data1$iraq, na.rm = TRUE)
	stem <- mean(data1$stem, na.rm = TRUE)
	wage <- mean(data1$wage, na.rm = TRUE)
	
## 2x2, true x party effects over party unity ##

	range <- seq(0.75, 1.00, 0.005)
	
	yesA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkA <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkB <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkC <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))

	yesD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	noD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))
	dkD <- matrix(rep(NA, 5000*length(range)), nrow = 5000, ncol = length(range))


## full model pics ##

	for(i in 1:length(range)){

	yes <- data$bbt1 * 0 +
		data$bbt3 * 1 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 0 * range[i] + 
		data$bbt15 * 0 * 0 + 
		data$bbt17 * 1 * range[i] + 
		data$bbt19 * 1 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 0 * polInt + 
		data$bbt31 * 1 * polInt + 
		data$bbt33 * 1 * polInt + 
		data$bbt35 * 0 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 1 * range[i] * polInt + 
		data$bbt41 * 1 * 0 * polInt + 
		data$bbt43 * 1 * range[i] * polInt + 
		data$bbt45 * 1 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 0 * 0.5 + 
		data$bbt55 * 1 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
	no <- data$bbt2 * 0 + 
		data$bbt4 * 1 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 0 * range[i] + 
		data$bbt16 * 0 * 0 + 
		data$bbt18 * 1 * range[i] + 
		data$bbt20 * 1 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 0 * polInt + 
		data$bbt32 * 1 * polInt + 
		data$bbt34 * 1 * polInt + 
		data$bbt36 * 0 * range[i] * polInt + 
		data$bbt38 * 0 * 0 * polInt + 
		data$bbt40 * 1 * range[i] * polInt + 
		data$bbt42 * 1 * 0 * polInt + 
		data$bbt44 * 1 * range[i] * polInt + 
		data$bbt46 * 1 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 0 * 0.5 + 
		data$bbt56 * 1 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84
		
	yesA[, i] <- yes
	noA[, i] <- no
	
	denom <- 1 + exp(yesA[, i]) + exp(noA[, i])
	yesA[, i] <- exp(yesA[, i]) / denom
	noA[, i] <- exp(noA[, i]) / denom
	dkA[, i] <- 1 - yesA[, i] - noA[, i]

	##########***********##########
	
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 1 * polInt + 
		data$bbt31 * 0 * polInt + 
		data$bbt33 * 1 * polInt + 
		data$bbt35 * 1 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 0 * range[i] * polInt + 
		data$bbt41 * 0 * 0 * polInt + 
		data$bbt43 * 1 * range[i] * polInt + 
		data$bbt45 * 1 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 1 * polInt + 
		data$bbt32 * 0 * polInt + 
		data$bbt34 * 1 * polInt + 
		data$bbt36 * 1 * range[i] * polInt + 
		data$bbt38 * 1 * 0 * polInt + 
		data$bbt40 * 0 * range[i] * polInt + 
		data$bbt42 * 0 * 0 * polInt + 
		data$bbt44 * 1 * range[i] * polInt + 
		data$bbt46 * 1 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84
		
	yesB[, i] <- yes
	noB[, i] <- no
	
	denom <- 1 + exp(yesB[, i]) + exp(noB[, i])
	yesB[, i] <- exp(yesB[, i]) / denom
	noB[, i] <- exp(noB[, i]) / denom
	dkB[, i] <- 1 - yesB[, i] - noB[, i]

	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 1 * polInt + 
		data$bbt31 * 0 * polInt + 
		data$bbt33 * 0 * polInt + 
		data$bbt35 * 1 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 0 * range[i] * polInt + 
		data$bbt41 * 0 * 0 * polInt + 
		data$bbt43 * 0 * range[i] * polInt + 
		data$bbt45 * 0 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 1 * polInt + 
		data$bbt32 * 0 * polInt + 
		data$bbt34 * 0 * polInt + 
		data$bbt36 * 1 * range[i] * polInt + 
		data$bbt38 * 1 * 0 * polInt + 
		data$bbt40 * 0 * range[i] * polInt + 
		data$bbt42 * 0 * 0 * polInt + 
		data$bbt44 * 0 * range[i] * polInt + 
		data$bbt46 * 0 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84

		
	yesC[, i] <- yes
	noC[, i] <- no
	
	denom <- 1 + exp(yesC[, i]) + exp(noC[, i])
	yesC[, i] <- exp(yesC[, i]) / denom
	noC[, i] <- exp(noC[, i]) / denom
	dkC[, i] <- 1 - yesC[, i] - noC[, i]

	##########***********##########
	
	
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * polInt + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * polInt + 
		data$bbt27 * 0 * polInt + 
		data$bbt29 * 1 * polInt + 
		data$bbt31 * 0 * polInt + 
		data$bbt33 * 1 * polInt + 
		data$bbt35 * 1 * range[i] * polInt + 
		data$bbt37 * 1 * 0 * polInt + 
		data$bbt39 * 0 * range[i] * polInt + 
		data$bbt41 * 0 * 0 * polInt + 
		data$bbt43 * 1 * range[i] * polInt + 
		data$bbt45 * 1 * 0 * polInt + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * polInt + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * polInt + 
		data$bbt28 * 0 * polInt + 
		data$bbt30 * 1 * polInt + 
		data$bbt32 * 0 * polInt + 
		data$bbt34 * 1 * polInt + 
		data$bbt36 * 1 * range[i] * polInt + 
		data$bbt38 * 1 * 0 * polInt + 
		data$bbt40 * 0 * range[i] * polInt + 
		data$bbt42 * 0 * 0 * polInt + 
		data$bbt44 * 1 * range[i] * polInt + 
		data$bbt46 * 1 * 0 * polInt + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84

	yesD[, i] <- yes
	noD[, i] <- no
	
	denom <- 1 + exp(yesD[, i]) + exp(noD[, i])
	yesD[, i] <- exp(yesD[, i]) / denom
	noD[, i] <- exp(noD[, i]) / denom
	dkD[, i] <- 1 - yesD[, i] - noD[, i]
		
	}

pdf('a4.3.pdf', width = 10, height = 12)
par(mfrow = c(2, 1))

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Inference Probability',
		main = 'Holding True Vote at Yea and\nChanging Party Line from Nay to Yea',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesD[, i]-yesC[, i], 0.975)
		hiN[i] <- quantile(noD[, i]-noC[, i], 0.975)
		hiD[i] <- quantile(dkD[, i]-dkC[, i], 0.975)
		loY[i] <- quantile(yesD[, i]-yesC[, i], 0.025)
		loN[i] <- quantile(noD[, i]-noC[, i], 0.025)
		loD[i] <- quantile(dkD[, i]-dkC[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesD-yesC), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noD-noC), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkD-dkC), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkD-dkC)[51], 'DK ')
	text(1.01, colMeans(noD-noC)[51], 'Nay')
	text(1.01, colMeans(yesD-yesC)[51], 'Yea')
	text(0.74, colMeans(dkD-dkC)[1], 'DK ')
	text(0.74, colMeans(noD-noC)[1], 'Nay')
	text(0.74, colMeans(yesD-yesC)[1], 'Yea')

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Inference Probability',
		main = 'Holding Party Line at Yea and\nChanging True Vote from Nay to Yea',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesB[, i]-yesA[, i], 0.975)
		hiN[i] <- quantile(noB[, i]-noA[, i], 0.975)
		hiD[i] <- quantile(dkB[, i]-dkA[, i], 0.975)
		loY[i] <- quantile(yesB[, i]-yesA[, i], 0.025)
		loN[i] <- quantile(noB[, i]-noA[, i], 0.025)
		loD[i] <- quantile(dkB[, i]-dkA[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesB-yesA), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noB-noA), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkB-dkA), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkB-dkA)[51], 'DK ')
	text(1.01, colMeans(noB-noA)[51], 'Nay')
	text(1.01, colMeans(yesB-yesA)[51], 'Yea')
	text(0.74, colMeans(dkB-dkA)[1], 'DK ')
	text(0.74, colMeans(noB-noA)[1], 'Nay')
	text(0.74, colMeans(yesB-yesA)[1], 'Yea')

dev.off()

## Interest Effects ##
## this will change interest from low to high ##
## for the case of a party line matching the true vote ##
## as well as for a mistmatch ##

	for(i in 1:length(range)){

	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 1 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * 1 + 
		data$bbt27 * 0 * 1 + 
		data$bbt29 * 1 * 1 + 
		data$bbt31 * 0 * 1 + 
		data$bbt33 * 1 * 1 + 
		data$bbt35 * 1 * range[i] * 1 + 
		data$bbt37 * 1 * 0 * 1 + 
		data$bbt39 * 0 * range[i] * 1 + 
		data$bbt41 * 0 * 0 * 1 + 
		data$bbt43 * 1 * range[i] * 1 + 
		data$bbt45 * 1 * 0 * 1 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 1 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * 1 + 
		data$bbt28 * 0 * 1 + 
		data$bbt30 * 1 * 1 + 
		data$bbt32 * 0 * 1 + 
		data$bbt34 * 1 * 1 + 
		data$bbt36 * 1 * range[i] * 1 + 
		data$bbt38 * 1 * 0 * 1 + 
		data$bbt40 * 0 * range[i] * 1 + 
		data$bbt42 * 0 * 0 * 1 + 
		data$bbt44 * 1 * range[i] * 1 + 
		data$bbt46 * 1 * 0 * 1 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84
		
	yesA[, i] <- yes
	noA[, i] <- no
	
	denom <- 1 + exp(yesA[, i]) + exp(noA[, i])
	yesA[, i] <- exp(yesA[, i]) / denom
	noA[, i] <- exp(noA[, i]) / denom
	dkA[, i] <- 1 - yesA[, i] - noA[, i]

	##########***********##########
	
	yes <- data$bbt1 * 1 +
		data$bbt3 * 0 + 
		data$bbt5 * 1 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 3 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 1 * range[i] + 
		data$bbt23 * 1 * 0 + 
		data$bbt25 * range[i] * 3 + 
		data$bbt27 * 0 * 3 + 
		data$bbt29 * 1 * 3 + 
		data$bbt31 * 0 * 3 + 
		data$bbt33 * 1 * 3 + 
		data$bbt35 * 1 * range[i] * 3 + 
		data$bbt37 * 1 * 0 * 3 + 
		data$bbt39 * 0 * range[i] * 3 + 
		data$bbt41 * 0 * 0 * 3 + 
		data$bbt43 * 1 * range[i] * 3 + 
		data$bbt45 * 1 * 0 * 3 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 1 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 1 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 3 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 1 * range[i] + 
		data$bbt24 * 1 * 0 + 
		data$bbt26 * range[i] * 3 + 
		data$bbt28 * 0 * 3 + 
		data$bbt30 * 1 * 3 + 
		data$bbt32 * 0 * 3 + 
		data$bbt34 * 1 * 3 + 
		data$bbt36 * 1 * range[i] * 3 + 
		data$bbt38 * 1 * 0 * 3 + 
		data$bbt40 * 0 * range[i] * 3 + 
		data$bbt42 * 0 * 0 * 3 + 
		data$bbt44 * 1 * range[i] * 3 + 
		data$bbt46 * 1 * 0 * 3 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 1 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt66 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84
		
	yesB[, i] <- yes
	noB[, i] <- no
	
	denom <- 1 + exp(yesB[, i]) + exp(noB[, i])
	yesB[, i] <- exp(yesB[, i]) / denom
	noB[, i] <- exp(noB[, i]) / denom
	dkB[, i] <- 1 - yesB[, i] - noB[, i]
	
		
	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 1 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * 1 + 
		data$bbt27 * 0 * 1 + 
		data$bbt29 * 1 * 1 + 
		data$bbt31 * 0 * 1 + 
		data$bbt33 * 0 * 1 + 
		data$bbt35 * 1 * range[i] * 1 + 
		data$bbt37 * 1 * 0 * 1 + 
		data$bbt39 * 0 * range[i] * 1 + 
		data$bbt41 * 0 * 0 * 1 + 
		data$bbt43 * 0 * range[i] * 1 + 
		data$bbt45 * 0 * 0 * 1 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 1 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * 1 + 
		data$bbt28 * 0 * 1 + 
		data$bbt30 * 1 * 1 + 
		data$bbt32 * 0 * 1 + 
		data$bbt34 * 0 * 1 + 
		data$bbt36 * 1 * range[i] * 1 + 
		data$bbt38 * 1 * 0 * 1 + 
		data$bbt40 * 0 * range[i] * 1 + 
		data$bbt42 * 0 * 0 * 1 + 
		data$bbt44 * 0 * range[i] * 1 + 
		data$bbt46 * 0 * 0 * 1 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84
	
	yesC[, i] <- yes
	noC[, i] <- no
	
	denom <- 1 + exp(yesC[, i]) + exp(noC[, i])
	yesC[, i] <- exp(yesC[, i]) / denom
	noC[, i] <- exp(noC[, i]) / denom
	dkC[, i] <- 1 - yesC[, i] - noC[, i]
		
		
	yes <- data$bbt1 * 1 + 
		data$bbt3 * 0 + 
		data$bbt5 * 0 + 
		data$bbt7 * range[i] + 
		data$bbt9 * 0 + 
		data$bbt11 * 3 + 
		data$bbt13 * 1 * range[i] + 
		data$bbt15 * 1 * 0 + 
		data$bbt17 * 0 * range[i] + 
		data$bbt19 * 0 * 0 + 
		data$bbt21 * 0 * range[i] + 
		data$bbt23 * 0 * 0 + 
		data$bbt25 * range[i] * 3 + 
		data$bbt27 * 0 * 3 + 
		data$bbt29 * 1 * 3 + 
		data$bbt31 * 0 * 3 + 
		data$bbt33 * 0 * 3 + 
		data$bbt35 * 1 * range[i] * 3 + 
		data$bbt37 * 1 * 0 * 3 + 
		data$bbt39 * 0 * range[i] * 3 + 
		data$bbt41 * 0 * 0 * 3 + 
		data$bbt43 * 0 * range[i] * 3 + 
		data$bbt45 * 0 * 0 * 3 + 
		data$bbt47 * 0.5 + 
		data$bbt49 * 0.5 + 
		data$bbt51 * 0.5 + 
		data$bbt53 * 1 * 0.5 + 
		data$bbt55 * 0 * 0.5 + 
		data$bbt57 * 0 * 0.5 + 
		data$bbt59 * 0.5 * 0.5 + 
		data$bbt61 * 0.5 * 0.5 + 
		data$bbt63 * female + 
		data$bbt65 * race + 
		data$bbt67 * income + 
		data$bbt69 * edu + 
		data$bbt71 * abort + 
		data$bbt73 * gains + 
		data$bbt75 * imm + 
		data$bbt77 * iraq + 
		data$bbt79 * stem + 
		data$bbt81 * wage + 
		data$bbt83
		
		
	no <- data$bbt2 * 1 + 
		data$bbt4 * 0 + 
		data$bbt6 * 0 + 
		data$bbt8 * range[i] + 
		data$bbt10 * 0 + 
		data$bbt12 * 3 + 
		data$bbt14 * 1 * range[i] + 
		data$bbt16 * 1 * 0 + 
		data$bbt18 * 0 * range[i] + 
		data$bbt20 * 0 * 0 + 
		data$bbt22 * 0 * range[i] + 
		data$bbt24 * 0 * 0 + 
		data$bbt26 * range[i] * 3 + 
		data$bbt28 * 0 * 3 + 
		data$bbt30 * 1 * 3 + 
		data$bbt32 * 0 * 3 + 
		data$bbt34 * 0 * 3 + 
		data$bbt36 * 1 * range[i] * 3 + 
		data$bbt38 * 1 * 0 * 3 + 
		data$bbt40 * 0 * range[i] * 3 + 
		data$bbt42 * 0 * 0 * 3 + 
		data$bbt44 * 0 * range[i] * 3 + 
		data$bbt46 * 0 * 0 * 3 + 
		data$bbt48 * 0.5 + 
		data$bbt50 * 0.5 + 
		data$bbt52 * 0.5 + 
		data$bbt54 * 1 * 0.5 + 
		data$bbt56 * 0 * 0.5 + 
		data$bbt58 * 0 * 0.5 + 
		data$bbt60 * 0.5 * 0.5 + 
		data$bbt62 * 0.5 * 0.5 + 
		data$bbt64 * female + 
		data$bbt64 * race + 
		data$bbt68 * income + 
		data$bbt70 * edu + 
		data$bbt72 * abort + 
		data$bbt74 * gains + 
		data$bbt76 * imm + 
		data$bbt78 * iraq + 
		data$bbt80 * stem + 
		data$bbt82 * wage + 
		data$bbt84
	
	yesD[, i] <- yes
	noD[, i] <- no
	
	denom <- 1 + exp(yesD[, i]) + exp(noD[, i])
	yesD[, i] <- exp(yesD[, i]) / denom
	noD[, i] <- exp(noD[, i]) / denom
	dkD[, i] <- 1 - yesD[, i] - noD[, i]
		
	}
	
## create the plots ##

pdf('a4.4.pdf', width = 10, height = 12)
par(mfrow = c(2, 1))

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Inference Probability',
		main = '',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesB[, i]-yesA[, i], 0.975)
		hiN[i] <- quantile(noB[, i]-noA[, i], 0.975)
		hiD[i] <- quantile(dkB[, i]-dkA[, i], 0.975)
		loY[i] <- quantile(yesB[, i]-yesA[, i], 0.025)
		loN[i] <- quantile(noB[, i]-noA[, i], 0.025)
		loD[i] <- quantile(dkB[, i]-dkA[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesB-yesA), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noB-noA), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkB-dkA), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkB-dkA)[51], 'DK ')
	text(1.01, colMeans(noB-noA)[51], 'Nay')
	text(1.01, colMeans(yesB-yesA)[51], 'Yea')
	text(0.74, colMeans(dkB-dkA)[1]+0.025, 'DK ')
	text(0.74, colMeans(noB-noA)[1]-0.025, 'Nay')
	text(0.74, colMeans(yesB-yesA)[1], 'Yea')
	mtext('Effect of Political Interest', line = 2.5, cex = 1.5, font = 2)
	mtext('Holding True Vote at Yea and Party Line at Yea\nChanging Political Interest from Low to High', line = 0.25)

	plot(1, 1, type = 'n',
		xlab = 'Past Unity',
		ylab = 'Change in Predicted Inference Probability',
		main = '',
		xlim = c(0.73, 1.02),
		ylim = c(-.6, .6),
		axes = FALSE)
	axis(1, at = seq(0.75, 1, 0.05))
	axis(2)
	polygon(c(0.75, 0.75, 1, 1), c(-0.6, 0.6, 0.6, -0.6), col = rgb(0, 0, 0, 0.1), border = NA)
	for(i in -6:6){
		lines(c(0.75, 1), c(i/10, i/10), col = 'white', lwd = 0.75)
	}
	for(i in 1:9){
		lines(c(0.75+(0.025*i), 0.75+(0.025*i)), c(-0.6, 0.6), col = 'white', lwd = 0.75)
	}
	hiY <- c()
	hiN <- c()
	hiD <- c()
	loY <- c()
	loN <- c()
	loD <- c()
	for(i in 1:length(range)){
		hiY[i] <- quantile(yesD[, i]-yesC[, i], 0.975)
		hiN[i] <- quantile(noD[, i]-noC[, i], 0.975)
		hiD[i] <- quantile(dkD[, i]-dkC[, i], 0.975)
		loY[i] <- quantile(yesD[, i]-yesC[, i], 0.025)
		loN[i] <- quantile(noD[, i]-noC[, i], 0.025)
		loD[i] <- quantile(dkD[, i]-dkC[, i], 0.025)
	}
	lines(c(0.75, 1), c(0, 0), col = rgb(0, 1, 0, 0.5), lwd = 2)
	polygon(c(range, rev(range)), c(hiD, rev(loD)), col = rgb(0, 0, 0, .1), border = NA)
	polygon(c(range, rev(range)), c(hiY, rev(loY)), col = rgb(0, 0, 1, .1), border = NA)
	polygon(c(range, rev(range)), c(hiN, rev(loN)), col = rgb(1, 0, 0, .1), border = NA)
	lines(range, colMeans(yesD-yesC), col = rgb(0, 0, 1, 0.25), lwd = 2)
	lines(range, colMeans(noD-noC), col = rgb(1, 0, 0, 0.25), lwd = 2)
	lines(range, colMeans(dkD-dkC), col = rgb(0, 0, 0, 0.25), lwd = 2)
	text(1.01, colMeans(dkD-dkC)[51], 'DK ')
	text(1.01, colMeans(noD-noC)[51], 'Nay')
	text(1.01, colMeans(yesD-yesC)[51], 'Yea')
	text(0.74, colMeans(dkD-dkC)[1], 'DK ')
	text(0.74, colMeans(noD-noC)[1], 'Nay')
	text(0.74, colMeans(yesD-yesC)[1], 'Yea')
	mtext('Effect of Political Interest', line = 2.5, cex = 1.5, font = 2)
	mtext('Holding True Vote at Yea and Party Line at Nay\nChanging Political Interest from Low to High', line = 0.25)
dev.off()

